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surface  normal  direction 

e 

dilatation  tensor 
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thermal  conductivity 
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Clauser’s  constant,  0.0168 
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molecular  and  turbulent  viscosities 
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molecular  and  turbulent  dynamic  viscosities 
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turbulence  model  viscosity  parameter 

i 

phase  angle 

n 
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a 
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T 

pseudo- time 

T 

viscous  stress  tensor 
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blade  row  angular  velocity 
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Q 
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Subscripts 
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i 

inviscid 
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max 

maximum 
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Chapter  1 


INTRODUCTION 

Flow  unsteadiness  is  the  dominant  aspect  of  flow  through  turbomachinery  compo¬ 
nents.  Examples  of  this  unsteadiness  include  turbulence  effects,  blade  row  interac¬ 
tions,  and  rotating  stall  and  surge.  Of  these,  blade  row  interaction  effects  have  the 
greatest  impact  on  the  performance  of  both  compressors  and  turbines.  Standard 
design  practice  ignores  much  of  this  unsteadiness,  using  steady-state  analyses  to 
predict  the  time-averaged  mean  behavior  of  the  flow.  As  design  goals  become  more 
aggressive,  it  becomes  increasingly  difficult  to  improve  performance  using  methods 
based  on  this  time-averaged  flow  assumption.  Therefore,  a  more  sophisticated  treat¬ 
ment  of  the  unsteady  nature  of  the  flow  becomes  vital.  Both  computational  and 
experimental  resources  need  to  be  used  to  characterize,  understand,  and  model  these 
effects  to  improve  the  capability  and  accuracy  of  current  design  tools.  The  ultimate 
use  for  a  computer  code  such  as  the  one  presented  in  this  thesis  is  to  provide  an 
insight  into  the  physics  of  the  interaction  between  moving  blade  rows. 

1.1  Background 

Current  turbomachinery  design  methods  use  the  S1-S2  streamsurface  method  first 
described  by  Wu  [80,  81],  with  two  intersecting  planes  used  to  describe  the  three- 
dimensional  space  within  a  turbomachinery  blade  passage  (see  Fig.  1.1).  Typically, 
an  entire  multi-stage  turbomachinery  compressor  or  turbine  is  designed  along  the 
meridional  r  —  z  plane  (Si  surface)  with  blade-to-blade  streamsurface  (S2  surface) 
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SI  surface 


Figure  1.1:  S1-S2  streamsurface  model  of  Wu. 

analyses  used  to  improve  the  performance  of  the  individual  cascade  sections.  In  order 
to  perform  the  meridional  analysis,  each  blade  row  is  designed  as  an  isolated  cascade. 
The  effect  of  the  up-  and  downstream  rows  is  only  included  as  a  circumferentially- 
averaged  “force  field.”  The  obvious  disadvantage  of  this  approach  is  that  it  removes 
the  impact  of  any  of  the  natural  flow  unsteadiness  from  the  design  methodology. 

An  example  of  the  impact  of  this  unsteadiness  may  be  found  in  transonic 
compressors.  In  this  case,  the  position  of  the  inlet  shock  pattern  determines  the  mass 
flow  through  the  blade  row  over  a  large  portion  of  the  operating  line.  The  transonic 
compressor  characteristic  is  typically  very  steep;  small  variations  in  mass  flow  may 
result  in  large  changes  in  the  operating  pressure  ratio.  Therefore,  unsteadiness  in 
the  shock  position  may  cause  the  blade  row  to  overflow  resulting  in  a  significant  drop 
in  the  pressure  ratio  from  the  design  point,  as  illustrated  in  Fig.  1.2.  This  effect  has 
been  observed  experimentally  by  Fleeter  et  al.  [25]  during  the  test  of  a  five-stage 
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compressor  which  had  fully  supersonic  second  and  third  stages.  They  observed  an 
approximate  5%  overflow  in  the  second  stage  from  the  design  value,  with  a  large 
drop  in  machine  pressure  ratio. 

Similar  effects  exist  in  transonic  turbine  designs.  The  two  main  concerns  are 
the  impingement  of  the  vane  trailing  edge  shock  on  the  downstream  rotor  and  the 
effect  of  unsteadiness  on  blade  heat  transfer.  Both  shock  impingement,  and  the 
accumulation  of  low  momentum  wake  fluid  along  the  rotor  pressure  surface,  can 
affect  the  heat  transfer  and  reduce  the  effectiveness  of  the  cooling  flows  added  to 
protect  the  blade  surfaces  from  the  high  temperatures  of  the  turbine  environment. 

The  purpose  of  this  thesis  is  to  present  a  numerical  technique  for  modeling  the 
unsteady  interaction  within  modern  turbomachinery  components  in  a  computation¬ 
ally  cost  effective  manner.  This  work  will  discuss  the  development  of  the  computer 
program,  its  implementation,  and  will  present  a  demonstration  of  its  performance 
on  a  modern  turbine  design. 

1.2  Blade  Row  Interactions 

Computationally,  the  blade  row  interaction  problem  has  been  studied  using  two 
major  approaches.  The  first  is  to  analyze  only  a  single  blade  row  while  imposing 
inlet  boundary  conditions  containing  a  representation  of  the  blade  wakes  from  an 
upstream  blade  row.  Mitchell  [49]  and  Hodson  [31]  used  this  approach  to  extend 
the  inviscid  time-marching  method  of  Denton  [24]  to  investigate  the  blade  surface 
static  pressure,  and  velocity  amplitude  and  phase  relationships,  at  the  wake  passing 
frequency.  Giles  [26]  solved  the  unsteady  Euler  equations  using  a  lagged  periodicity 
boundary  condition  technique,  which  inclines  the  computational  plane  in  time  to 
allow  calculation  of  arbitrary  rotor /stator  blade  count  ratios.  Scott  and  Hankey  [66] 
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Figure  1.2:  Blade  row  interaction  effect  on  transonic  compressor  design  point. 
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and  the  present  author  [59]  solved  the  Navier-Stokes  equations  in  efforts  to  calculate 
the  unsteady  mass  flow  swallowing  capacity  of  a  compressor  rotor  cascade. 

The  second  approach  is  to  model  both  the  rotor  and  stator  blade  rows  as  a 
fully  coupled  system,  with  the  blade  wakes  generated  by  an  upstream  row  impinging 
on  the  downstream  blades.  This  approach  has  the  added  advantage  that  both  shock 
impingement  and  the  potential  interaction  between  blades  is  modeled,  as  well  as 
any  changes  in  the  upstream  wake  due  to  different  stage  loading  conditions.  The 
first  two-dimensional  calculations  were  done  by  Rai  [52]  using  a  third-order  accurate 
implicit  scheme  with  patched  and  overlaid  grids  to  solve  the  thin-layer  Navier-Stokes 
equations  in  a  low-speed  turbine  stage.  This  method  was  later  extended  to  compute 
a  low-speed,  multi-stage  compressor  by  Gundy-Burlet  et  al.  [30]  and  a  multi-stage 
turbine  by  Lin  and  Yang  [42],  Other  two-dimensional  codes  have  included  the  explicit 
codes  of  Jorgenson  and  Chima  [36],  Lewis  et  al.  [41],  Arnone,  et  al.  [4,  5],  and 
the  present  author  [60].  Three-dimensional  calculations  were  first  presented  by  Rai 
[53,  54]  for  a  low-speed  turbine  which  included  the  effect  of  tip  leakage  over  the 
rotor.  Other  prominent  three-dimensional  codes  include  the  set  of  inviscid  and 
viscous  codes  of  Rao  and  Delaney  [56]  and  Rao  et  al.  [57,  58].  These  methods  all  use 
multiple  structured  meshes,  generally  with  patched  or  Chimera-type  grid  structures. 
Dawes  [23]  and  Mathur  et  al.  [44,  45]  have  recently  developed  unstructured  codes  to 
compute  rotor-stator  interactions,  but  neither  of  these  methods  includes  adaptation 
of  the  grid  to  the  unsteady  flow  field. 

1.3  Numerical  Technique 

Generally,  methods  used  for  the  blade-row  interaction  problem  grid  each  blade  row 
separately  and  then  use  a  sliding  interface  to  transfer  information  between  the  dif¬ 
ferent  blade  rows.  In  the  case  of  modern  high-speed  compressor  and  turbine  designs, 
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the  flow  fields  contain  unsteady  shock  and  wake  patterns  which  cross  this  interface 
between  grids.  A  large  change  in  grid  density  across  the  interface  will  induce  a 
significant  dissipation  of  both  wakes  and  shocks  causing  damping  of  the  interaction 
effects.  To  adequately  resolve  the  complex  shock  and  wake  structures  using  struc¬ 
tured  grids  requires  a  considerably  larger  number  of  points  than  is  typically  used, 
most  of  them  wasted  in  other,  more  benign,  regions  of  the  flow. 

The  use  of  unstructured  grids  is  an  alternative  approach  that  removes  many  of 
the  restrictions  found  in  the  structured  grid  approach.  Virtually  any  grid  point  dis¬ 
tribution  can  be  used,  allowing  grid  points  to  be  added  only  where  there  is  a  need  for 
increased  resolution.  Coupled  with  a  solution-adaptive  algorithm,  unstructured  grids 
allow  the  computational  grid  to  evolve  with  the  solution,  reducing  both  computer 
time  and  storage  requirements  over  those  for  a  structured  grid  of  adequate  resolu¬ 
tion.  Recent  literature  has  shown  a  number  of  examples  where  unstructured  grids 
have  been  used  for  steady  [61,  69]  as  well  as  unsteady  [29,  44,  45]  turbomachinery 
flows.  There  are  few  examples  in  the  literature,  however,  where  solution-adaptivity 
has  been  used  in  the  computation  of  unsteady  rotor-stator  flows. 

Another  factor  in  the  choice  of  a  computational  method  is  the  selection  of  an 
explicit  or  implicit  algorithm.  Implicit  solvers  offer  the  benefit  of  larger  time  steps, 
which  are  particularly  attractive  for  the  highly  stretched  grids  needed  to  adequately 
resolve  viscous  layers.  For  unsteady  calculations,  where  small  time  steps  are  often 
required  for  adequate  temporal  resolution,  the  computational  cost  of  an  implicit 
solver  used  throughout  the  entire  domain  may  overshadow  the  solver’s  ability  to 
operate  with  larger  time  steps.  Further,  implicit  solvers  for  unstructured  meshes 
are  computationally  quite  expensive.  Jameson  [32]  and  Weiss  and  Smith  [76]  have 
introduced  similar  techniques  which  use  an  explicit  time  marching  method  as  the 
driver  for  a  fully  implicit  scheme.  These  methods,  winch  are  adopted  in  this  thesis, 
provide  a  considerable  improvement  over  a  purely  explicit  implementation.  The 
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technique  of  Jameson  has  been  applied  to  unsteady  airfoil  flows  [32]  and  recently  by 
Arnone  for  blade  row  interactions  in  turbomachinery  [4,  5]. 

1.4  Turbulence  Models 

One  of  the  requirements  for  the  turbulence  model  used  in  this  effort  was  that  it 
should  contain  some  form  of  the  field  equations  to  allow  the  propagation  of  the 
eddy  viscosity  from  an  upstream  to  a  downstream  blade  row.  A  number  of  models 
are  currently  available,  both  in  one-  and  two-equation  form.  Two-equation  model 
choices  include  those  of  Jones  and  Launder  [35],  Coakley  [16],  Wilcox  [79],  and 
Menter  [48].  Of  these,  the  blended  k  —  e/k  —  u)  of  Menter  has  perhaps  the  strongest 
record  for  the  prediction  of  adverse  pressure  gradient  flows.  However,  the  additional 
computational  overhead  of  a  sixth  field  equation  led  to  a  consideration  of  the  one- 
equation  models  of  Johnson  and  King  [34],  Baldwin  and  Barth  [8],  and  Spalart  and 
Allmaras  [71].  The  Spalart  and  Allmaras  model  was  chosen  for  the  current  effort; 
it  has  the  advantage  of  being  a  local  model,  easily  implemented  in  an  unstructured 
grid  environment,  and  is  particularly  well  suited  to  predictions  with  coarser  grids. 

1.5  Adaptation  Methods 

Solution-adaptation  techniques  are  generally  classified  into  three  basic  categories. 
The  first,  and  one  that  may  be  applied  to  either  structured  or  unstructured  meshes, 
is  mesh  movement  in  response  to  flow  field  topology  [11,  20].  Original  mesh  connec¬ 
tivity  is  maintained,  with  points  moving  towards  or  away  from  regions  of  high  or 
low  gradients,  respectively.  The  obvious  disadvantage  of  this  approach  is  the  restric¬ 
tion  to  the  original  number  of  mesh  points  and  their  original  connectivity.  Complex 
flow  features  may  rob  the  remaining  mesh  of  resolution,  resulting  in  an  overall  grid 
distribution  which  is  unsatisfactory. 
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The  second  approach,  and  the  one  primarily  used  with  triangular  element 
meshes,  is  mesh  regeneration,  where  the  computational  mesh  is  completely  renewed 
during  an  adaptation  cycle  [50].  One  advantage  of  this  approach  is  that  no  change 
is  required  to  the  flow  solver  since  the  logical  mesh  structure  is  not  modified  after 
refinement.  Additionally,  re-meshing  gives  superior  control  over  variations  in  grid 
cell  size  and  allows  alignment  of  the  grid  with  flow  features.  The  primary  drawback 
of  this  method  is  that  the  computational  time  required  to  perform  the  re-meshing 
may  be  prohibitive,  particularly  when  the  mesh  must  be  redefined  repeatedly  as  in 
the  case  of  a  rotor-stator  interaction  type  problem. 

The  third  choice,  and  the  one  used  in  the  current  solver,  is  grid  enrichment 
through  local  cell  division  and  recombination,  the  technique  most  commonly  used 
with  quadrilateral  meshes.  Here  the  primary  advantage  is  the  speed  at  which  the 
grid  may  be  adapted,  of  prime  importance  in  the  computation  of  unsteady  flows. 
Special  coding  logic  is  required  to  account  for  the  midface  nodes  created  at  the 
interface  between  divided  and  undivided  cells.  A  series  of  rules  is  applied  during  the 
adaptation  to  insure  that  cells  divide  in  the  proper  sequence  to  inhibit  the  formation 
of  cells  with  multiple  midface  nodes  on  a  side.  Triangular  cell  unstructured  mesh 
methods  may  also  use  local  cell  refinement,  but  the  refinement  must  be  followed  by 
a  local  retriangulation  to  insure  a  smooth  mesh  [9].  It  is  unclear  from  the  literature 
whether  there  is  a  performance  advantage  between  quadrilateral  and  triangular  grid 
enrichment  techniques. 


1.6  Present  Research 

The  objective  of  the  current  study  is  the  development  of  an  unsteady  solution- 
adaptive  unstructured  grid  method  for  application  to  blade  row  interactions  in  tran¬ 
sonic  turbomachinery.  While  the  method  is  primarily  for  axial  cascades,  the  equa- 
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tions  are  cast  in  the  more  general  quasi-three-dimensional  formulation  to  include  the 
effects  of  radius  and  streamtube  thickness  variation  found  in  centrifugal  and  mixed 
flow  designs. 

The  algorithm  used  to  perform  the  time  integration  of  the  flow  equations  uses 
a  modification  of  the  implicit  reformulation  of  Jameson  [32],  allowing  time  accurate 
computations  and  the  use  of  convergence  acceleration  techniques.  The  effects  of 
turbulence  is  modeled  with  the  either  the  algebraic  model  of  Baldwin  and  Lomax 
[7]  or  the  one-equation  model  of  Spalart  and  Allmaras  [71]. 

Results  are  presented  for  a  variety  of  test  cases  used  to  separately  validate  the 
various  components  of  the  computer  code.  The  fully  unsteady  method  is  demon¬ 
strated,  both  with  and  without  grid  adaptation,  with  the  calculation  of  a  transonic 
turbine  design. 


1.7  Retrospective 

The  work  described  in  this  thesis  is  the  culmination  of  approximately  ten  years 
of  effort  in  the  area  of  modeling  of  unsteady  flows  in  turbomachinery.  This  section 
describes  a  number  of  the  directions  pursued  and  attempts  to  summarize  the  various 
methods  and  the  motivation  for  their  use. 

The  author’s  earliest  work  in  this  area  was  the  use  of  modified  version  of  the 
highly  successful  structured  code  developed  by  Shang  et  al.  [67,  68],  which  is  based 
on  the  explicit  MacCormack  predictor-corrector  method  [43].  This  code  was  used  by 
this  author  [59],  as  well  as  by  Scott  and  Hankey  [66],  to  predict  the  unsteady  mass 
flow  swallowing  capacity  of  a  transonic  fan  cascade  using  an  unsteady  inlet  wake 
boundary  condition.  Although  these  efforts  did  show  an  effect  of  the  upstream  wake 
parameters  on  the  time  averaged  mass  flow,  questions  remained  on  the  validity  of 
the  wake  boundary  condition.  Additionally,  the  potential  flow  interaction  between 


blade  rows  was  not  included  in  the  analysis.  These  considerations  led  to  a  further 
modification  of  the  structured  code  to  predict  a  fully  interactive  rotor-stator  flow 
field. 

While  this  code  was  successful  in  computing,  to  the  author’s  knowledge,  the 
first  prediction  of  rotor-stator  interaction  in  a  transonic  compressor  [60],  it  suffered 
from  a  number  of  deficiencies.  The  first  of  these  was  the  use  of  a  grid  with  inadequate 
density.  When  computing  an  isolated  airfoil  or  cascade  flow,  the  usual  practice  is  to 
stretch  the  grid  near  the  blade  surface  and  downstream  of  the  trailing  edge  in  order  to 
adequately  resolve  the  large  velocity  gradients  in  the  blade  boundary  layer  and  wake 
regions.  Typically,  the  remainder  of  the  flow  domain  is  relatively  sparsely  gridded 
because  the  gradients,  except  in  the  case  of  shocks,  are  considerably  more  mild. 
Correct  resolution  of  the  blade  wake  is  of  particular  importance  in  turbomachinery 
cascades  for  proper  prediction  of  cascade  losses  in  both  compressors  and  turbines, 
as  well  as  the  blade  heat  transfer  in  turbines.  Shock  resolution  has  also  received 
considerable  attention  in  external  airflow  calculations,  with  airfoil  grids  often  clus¬ 
tered  in  the  regions  of  expected  shocks,  in  an  effort  to  improve  the  resolution  in  the 
shock-boundary  layer  interaction  region.  In  turbomachinery,  shocks  play  a  consider¬ 
ably  larger  role:  setting  the  mass  flow  swallowing  capacity  of  transonic  compressors, 
and  impacting  transonic  turbine  heat  transfer  and  aeromechanical  response. 

Including  the  propagation  of  upstream  blade  wakes  into  a  downstream  row  is 
important  for  the  modeling  of  the  unsteady  flow  field.  To  accomplish  this  requires 
that  the  grid  resolution  is  maintained  between  up-  and  downstream  rows.  Since  the 
downstream  row  is  moving  relative  to  the  upstream,  the  circumferential  grid  resolu¬ 
tion  in  the  downstream  row  would  need  to  everywhere  match  the  finest  resolution  in 
the  upstream  row.  A  similar  requirement  exists  for  the  grid  clustering  needed  around 
the  shocks  that  propagate  between  blade  rows.  The  need  for  significantly  enhanced 
grid  resolution  obviously  poses  a  considerable  challenge  to  both  memory  and  com- 
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puting  time  and  was  the  prime  driver  behind  the  selection  of  a  solution-adaptive 
unstructured  grid  method. 

Another  issue  in  the  original  rotor-stator  code  was  the  use  of  an  algebraic 
turbulence  model,  the  industry  standard  at  the  time.  An  algebraic  model  computes 
the  turbulent  eddy  viscosity  by  searching  away  from  the  wall  surface  to  determine 
which  type  of  profile  model  to  apply.  This  search  is  often  computationally  expensive 
in  unstructured  grids,  which  do  not  have  the  natural  connectivity  between  grid  points 
contained  in  structured  grids.  More  important  however,  any  method  to  model  the 
propagation  of  the  the  turbulent  eddy  viscosity  from  one  blade  row  to  another  would 
be  woefully  inadequate  using  an  algebraic  model.  Therefore,  it  was  felt  that  any 
turbulence  model  used  for  rotor-stator  calculations  must  include  some  form  of  the 
field  equations  in  order  to  capture  the  effect  of  turbulent  eddy  viscosity  convection 
and  diffusion.  An  added  difficulty  with  the  algebraic  model  was  its  extremely  poor 
performance  in  the  regions  of  strong  pressure  gradient  often  found  in  modern,  highly 
loaded,  transonic  compressors.  Calculations  showed  an  unrealistically  large  region 
of  separation,  causing  a  large  wake  which  would  create  a  “subsonic  gap”  in  the  flow- 
entering  the  rotor.  The  unsteady  animation  of  this  flow  showed  a  section  of  the  shock 
wave  propagating  forward  inside  the  low-speed  flow  of  the  stator  wake,  triggering  an 
axial  oscillation  of  the  stator  separation  point.  These  computations  clearly  showed  a 
secondary  pulsation  of  the  flow  caused  by  the  shock  induced  unsteadiness  within  the 
wake.  While  quite  entertaining  to  watch,  there  is  little  doubt  that  the  inadequacies 
of  the  turbulence  model  were  largely  responsible  for  this  phenomenon.  An  ability 
to  properly  account  for  strong  pressure  gradients  is  therefore  also  included  in  the 
turbulence  model  selection  criteria. 
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1.8  Overview  of  Thesis 


In  Chapter  2,  the  Navier-Stokes  equations  are  derived  in  a  general  sense,  including  a 
description  of  all  assumptions  required  for  the  present  effort.  Using  the  more  detailed 
analysis  included  in  Appendix  A,  these  equations  are  presented  in  the  specific  quasi- 
3D  form,  which  includes  the  effects  of  radius  and  streamtube  thickness  variation. 
The  method  of  incorporating  these  distributions  is  then  examined,  followed  by  an 
explanation  of  the  nondimensionalization  of  the  flow  equations. 

Chapter  3  describes  the  finite  volume  solution  procedure  implemented  on  an 
unstructured  mesh  of  quadrilateral  cells.  The  method  used  for  the  evaluation  of  the 
inviscid  and  viscous  terms  is  presented,  with  a  discussion  of  the  boundary  conditions 
and  the  inclusion  of  artificial  dissipation.  The  dual  time  stepping  algorithm  and 
convergence  acceleration  techniques  are  described,  followed  by  an  examination  of  the 
numerical  properties  of  the  scheme.  The  chapter  is  completed  with  a  description  of 
the  implementation  of  the  computer  program  and  its  semi-structured  grid  approach. 

Chapter  4  presents  a  discussion  of  the  two  turbulence  models  available  in  the 
code.  The  simple  algebraic  model  of  Baldwin  and  Lomax  [7]  is  first  described, 
followed  by  the  one-equation  model  of  Spalart  and  Allmaras  [71].  A  discussion  of 
the  special  requirements  of  the  unstructured  grid  implementations  of  both  methods 
is  presented. 

In  Chapter  5  the  method  used  for  the  generation  and  adaptation  of  the  unstruc¬ 
tured  grid  is  described.  Cell  division  and  recombination  is  used  to  either  add  or 
remove  individual  cells  from  the  mesh  using  the  local  flow  gradients  to  drive  the 
addition  or  removal  of  mesh  points.  Adaptation  of  the  grid  to  both  strong  and  weak 
flow  features  is  described,  followed  by  an  examination  of  the  numerical  treatment  of 
the  midface  nodes. 


12 


Chapter  6  presents  the  results  of  the  numerical  experiments  performed  for  this 
effort.  Validation  cases  are  presented  which  exercise  various  portions  of  the  code 
independently,  followed  by  a  sample  transonic  turbine  application.  A  number  of 
turbine  cases  are  presented  to  highlight  the  consequences  of  the  various  parameter 
choices  available  in  the  computation  of  unsteady  flows. 

Finally,  Chapter  7  presents  the  major  conclusions  of  this  work.  This  chapter 
also  includes  a  discussion  of  the  limitations  of  this  effort  and  recommendations  for 
future  work. 

The  thesis  is  completed  with  three  appendices.  The  first  presents  a  full  deriva¬ 
tion  of  the  quasi-3D  form  of  the  flow  equations.  The  second  contains  a  detailed 
description  of  the  boundary  conditions  implemented  in  the  computer  program.  The 
interactive  grid  generator  used  to  construct  the  computational  grids  is  described  in 
the  final  appendix. 
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Chapter  2 


EQUATIONS  OF  MOTION 


In  this  chapter  the  general  Navier-Stokes  equations  are  presented,  including  the 
assumptions  that  are  required  for  their  derivation.  The  quasi-three-dimensional  or 
quasi-3D  formulation  of  these  equations  is  then  given  and  expressed  in  a  Reynolds 
averaged,  nondimensional  form. 

2.1  Navier-Stokes  Equations 

Consider  the  two-dimensional  representation  of  an  arbitrary  control  volume  V  as 

shown  in  Fig.  2.1.  This  volume  is  bounded  by  surface  Sv,  with  outward  normal 

vector  n.  If  this  bounding  surface  moves  with  velocity  Vs  and  the  fluid  velocity 

— ^ 

relative  to  the  surface  is  Vr,  then  the  absolute  velocity  of  the  fluid  is  given  by 
V  =  Vs  +  Vr. 

Conservation  of  mass  within  this  control  volume  requires  that  the  time  rate- 
of-change  of  mass  within  the  volume  is  balanced  by  the  net  flux  of  mass  across  the 
bounding  surfaces,  i.e.: 

it  fJf^dv  +  i^-MS  =  °-  w 

This  equation  is  based  on  the  assumption  that  the  fluid  is  continuous  and  that  all 
flow  variables  are  continuous  functions  in  space. 

Newton’s  second  law  requires  that  the  time  rate- of- change  of  momentum  within 
the  volume,  plus  the  momentum  flux  convected  through  the  boundary,  is  balanced 
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Sv 


Figure  2.1:  Control  volume  for  conservation  equations. 

by  the  sum  of  all  forces  acting  on  the  boundary  of  the  control  volume.  This  is 
expressed  as: 

7tJI IvP^rdV  +  /s  •  ”) dS  =  (2-2) 

Neglecting  body  force  terms,  the  forces  acting  on  the  boundary  are  expressed  as  a 
summation  of  pressure  and  shear  stress  terms: 

EF  =  l  (—p  +  r)  •  hdS,  (2.3) 

J  Sy 

where  p  is  the  local  static  pressure  and  r  is  the  shear  stress  tensor. 

The  first  law  of  thermodynamics  states  that  the  time  rate-of-change  of  the 
total  energy  E  within  the  control  volume,  plus  the  energy  flux  convected  through 
the  boundary,  is  balanced  by  the  sum  of  the  heat  flux  Q  and  work  W  done  by  the 
stresses  on  the  fluid.  This  relation  is  given  by: 

jJJJvEdV  +  jf  E(Vr  ■  n)dS  =  Q-W.  (2.4) 
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Where  Q  includes  heat  conduction  alone,  with  the  assumption  that  there  are  no 
heat  sources  internal  to  the  control  volume.  The  work  W  is  given  by: 

W  —  <f  (TV  —  pV)  ■  ndS.  (2.5) 

d  Sy 

Auxiliary  relationships  are  required  to  close  the  above  system  of  equations. 
The  first  is  the  assumption  that  the  working  fluid  is  a  thermally  perfect  gas: 

p  =  pRT,  (2.6) 

where  R  is  the  gas  constant.  Further,  the  gas  is  assumed  to  be  calorically  perfect 
(i.e.,  with  constant  specific  heats  cv  and  cp).  The  specific  internal  energy  of  the  gas 
is  then  defined  as: 

e  =  cvT.  (2.7) 

Using  the  relationship  between  the  static  pressure  p  and  the  total  energy  of  the  fluid 
E,  the  above  equations  can  be  combined  to  yield: 

E  ~  Pe  +  9  (um  +  vf)  (2-8) 

and 

P=(7-l)[*-§  (<  +  »?)]  .  (2.9) 

where  7  is  the  ratio  of  specific  heats,  and 

f7  'Um'l'm  "h  VqZq.  (^2.10) 
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The  viscosity  of  the  fluid  y  is  expressed  as  a  function  of  temperature  using 
Sutherland’s  Law: 

h  _  ( j£_  \  2  (-^0°  ~h  Tref)  (9  ii\ 

yoo~\Tco)  T  +  TreS )’  ’ 

where  the  reference  temperature  Trej  =  1 1 0.4 A'  for  air,  and  y  and  y03  are  the 
viscosities  at  temperatures  T  and  T^. 

Finally,  the  thermal  conductivity  of  the  working  fluid  k  is  specified  using  the 
Prandtl  number: 

Pr  =  (2.12) 

K 

A  value  of  Pr  —  0.72  is  normally  taken  for  air  at  standard  conditions. 

2.2  Quasi-3D  Formulation 

Until  this  point  in  the  derivation,  an  arbitrary  three-dimensional  control  volume  has 
been  used.  The  specific  control  volume  used  to  produce  the  equations  in  their  final 
quasi-3D  form  must  now  be  described.  The  configuration  of  interest  is  the  compu¬ 
tation  of  the  blade-to-blade  flow  in  axial  and  radial  turbomachinery  with  the  inclu¬ 
sion  of  radius  and  streamsheet  thickness  variations.  Therefore,  a  two-dimensional 
axisymmetric  coordinate  system  (m,  9)  is  used,  where  m  is  the  streamwise  direction 
and  6  is  the  cylindrical  coordinate  rotating  with  the  blade  row: 

dm2  =  dz 2  +  dr2,  (2-13) 

where  a  is  the  axial  direction,  and 

6  =  9'  -  Qt.  (2.14) 
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Figure  2.2:  Quasi-3D  streamsurface  through  a  typical  fan  rotor. 


Here  6l  is  fixed  in  space  and  fi  is  the  angular  velocity  of  the  blade  row.  Additionally, 
the  velocity  normal  to  the  streamtube  lamina  is  assumed  to  be  zero.  The  streamsur¬ 
face  for  this  coordinate  system  through  a  typical  fan  rotor  is  shown  in  Fig.  2.2.  The 
computational  domain  is  a  two-dimensional  blade-to-blade  surface  discretized  with 
a  grid  of  computational  elements  or  cells.  These  cells  have  a  thickness  and  radial 
position  which  is  accounted  for  in  the  quasi-3D  formulation. 

In  this  coordinate  system,  the  radius  r  and  the  streamtube  thickness  h  are 
taken  to  be  known  functions  of  m.  The  equations  presented  below  reduce  to  the 
standard  two-dimensional  Navier-Stokes  equations  in  conservative  form  when  r  and 
h  are  constants.  The  validation  cases  presented  later  are  computed  with  r  —  h  =  1, 
while  the  cascade  computations  are  performed  with  radius  and  streamtube  thickness 
variations  defined  from  experimental  or  design  data. 
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The  Navier-Stokes  equations  for  this  domain  may  be  expressed  in  an  integral 
formulation,  which  is  used  in  a  finite  volume  implementation  as  follows: 


—  f  f  „  U dmdO  +  I  c e,i  ( Fd9  -  Gdm )  = 

dt  J  J  area  J  faces 

[ [  ccll  Hdmde  +  <f  cell  ( Rd6  -  Sdm ).  (2.15) 

Here,  t  is  time,  U  is  the  dependent  variable  vector.  F  and  G  are  the  inviscid  flux 
vectors,  R  and  S  are  the  viscous  flux  vectors,  and  H  is  the  source  term  vector 
incorporating  the  additional  radius  and  streamtube  thickness  variation  terms.  The 
full  derivation  of  the  quasi-3D  equations  is  presented  in  Appendix  A.  The  resulting 
vector  components  are  given  by  (using  the  notation  of  Chima  [15]  ): 


,  F  —  rh 


PGa  +  P 
{pvmve)r 
Vm^F  T  p) 


G  =  h 


pweVm 


(pweve  +  p)r 
wg(E  +  p)  +  rOp 


H  =  rh 


R  =  rh 


S  =  h 


(2.16) 


^m^mm  ~i~  ^ 6  1  mO  Qm 


Vm  1  mO  +  VQ  /  QQ  (f& 
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The  source  term  is  given  by: 


H2  =  ( pvg  +p-Tee)—  +  (p-  Trr)^, 

r  a 


using  the  notation 


1  dr 
r  dm J 


hm  1  dh 
h  h  dm 


The  static  pressure  is  given  by: 


P  =  (  7-1) 


E 


(4  +  »?) 


(2.17) 


and  the  stress  and  heat  conduction  terms  are: 


Tmm  =  +  ^9, 


Tee  -  2p  (--Qjf  +  +  A0, 


,  dve  1  dv. 
TmO  T 


m  T'm  \ 

9m  '  r  d6  9  r  ’ 


2/mm  — — f  A0, 
h 


dT 


tym  —  ^  . 


dm5 


1  dT 

<18  ~  K-77, 
r  90 


(2,18) 


where  w$  =  v$  —  rO  is  the  relative  velocity,  A  =  —2^/3,  and  the  dilatation  is  given 
by: 

0  = 


Trr,  ,  ^'771 


9m  r  90  m  V  r  h 


2.3  Reynolds  Averaging 


Flows  in  modern  turbomachinery  components  are  highly  turbulent  due  to  the  high 
Reynolds  numbers  of  the  flows.  Turbulence  is  characterized  by  flow  unsteadiness  over 
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a  wide  range  of  time  and  length  scales,  the  resolution  of  which  is  quite  impractical  for 
actual  turbomachinery  geometries  using  current  computing  resources.  Instead,  the 
high  frequency  unsteadiness  of  turbulent  flow  is  modeled  by  assuming  that  all  flow 
variables  are  comprised  of  two  components.  The  first  is  a  value  which  is  Reynolds 
averaged  over  a  given  time  interval  and  the  second  component  which  fluctuates  about 
this  average  to  represent  the  turbulence.  The  Reynolds  averaging  (or  more  correctly, 
mass,  or  Favre  averaging  for  compressible  flows)  results  in  the  addition  of  apparent 
stress  and  heat  transfer  terms  to  the  momentum  and  energy  equations.  Generally, 
the  Bussinesq  approximation  is  used,  which  assumes  that  these  additional  terms 
are  proportional  to  the  mean  strain  rate  and  temperature  gradient,  respectively.  In 
essence,  this  modeling  takes  the  form  of  a  redefinition  of  the  viscosity  and  thermal 
conductivity: 

fteff  —  ^  d"  (2.19) 

where  //  and  k  represent  the  molecular  values  of  the  viscosity  and  thermal  con¬ 
ductivity  and  ht  and  kt  are  their  turbulent  counterparts.  The  turbulent  Prandtl 
number 

PrT  =  ^  (2.20) 

kt 

is  assumed  constant,  with  a  value  of  PrT  =  0.90.  The  final  form  of  the  Reynolds 
averaged  stress  and  heat  transfer  terms  is  given  by: 


1  mm 
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g(^  +  ^r) 


9  dVm  _  1  dve 
^  dm  r  d6 
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Jdv6 
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dvm 

J m  ’ 
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Table  2.1:  Summary  of  nondimensional  reference  quantities. 


Trr  =  |(^  +  n)  L  (2^  - 


dvm 

dm 


1  dvg 

r  09 


(2.21) 


and 


C[m  — 


JL  4. 

Pr  Prj)  dm 


Jf_  ,  ±t_\  ldT 

qe  \Pr  PrT)rde' 


(2.22) 


The  models  used  to  evaluate  the  turbulent  eddy  viscosity  fij  are  presented  in  further 
detail  in  Chapter  4. 


2.4  Nondimensionalization 

These  equations  are  made  dimensionless  by  selection  of  a  characteristic  length  L, 
inlet  freestream  static  density  p <*,,  inlet  velocity  u <*,,  inlet  static  temperature  Too, 
and  molecular  viscosity  Table  2.1  gives  a  summary  of  the  reference  quan¬ 

tities  used  for  nondimensionalization.  As  a  notational  convenience,  all  variables  will 
henceforth  be  assumed  to  represent  their  nondimensional  equivalents. 
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The  only  effects  of  the  nondimensionalization  are  in  the  viscous  stress  and  heat 


conduction  terms  which  now  become: 


1"  m  . 


_  2  1  (  dvm  1  dve 


2  I  f^ldvg  .  („rm  hm 


dvri 


™  =  3^  +  n)Rz  [PltF  +  V"‘{2'f~  h/  dm)’ 

1  (dve  ldvm  rm\ 

Tm‘  =  (/■  +  n)Tt  +  -w  -  my)  , 

2  (  .  1  /  (  hm  vm  \  dvm  1  dvg  \ 

Trr  =  3  ^  +  ^~Re  [Vm  \  I  ~  V)  ~  ~  r~W  )  ’ 


1 


_ _ J£_  ,  Pt  \  dl [ 

qm~  (7  -  1  )ReMl  \Pr  PrT)  dm. ’ 

.  _ _ 1 _ f  P  .  Pt  h  1  dT 

qe~  (7  -  l)ReMl  \Pr  PrT)  r  dd’ 

where  the  Reynolds  number  is: 


(2.23) 


Re  — 


P'U'qqP 

Pco 


(2.24) 


The  characteristic  length  L  is  typically  taken  as  a  unit  length  in  the  coordinate 
system.  The  Mach  and  Reynolds  numbers,  and  Re,  and  the  inlet  flow  angle  a ^ 
are  specified  through  the  boundary  conditions. 
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Chapter  3 


UNSTRUCTURED  GRID  ALGORITHM 


This  chapter  contains  a  description  of  the  numerical  method  used  in  the  current  effort 
and  details  of  its  implementation.  A  description  is  provided  for  the  discretization 
of  the  inviscid  and  viscous  terms,  and  the  boundary  conditions  employed.  Imple¬ 
mentation  details  of  the  artificial  dissipation  operator  and  the  numerical  integration 
scheme  are  given  next,  followed  by  an  examination  of  the  numerical  properties  of 
the  overall  unstructured  scheme.  The  chapter  concludes  with  a  discussion  of  the 
semi-structured  implementation  strategy  and  the  resulting  data  structures  used  in 
the  computer  code. 


3.1  Inviscid  Terms 

A  one-step  Lax  Wendroff-style  integration  scheme  is  used  to  discretize  the 
convective  terms  in  the  inviscid  portion  of  Eqn.  2.15  on  a  stationary  grid  of  quadri¬ 
lateral  cells.  Considering  a  single  cell,  as  in  Fig.  3.1,  the  equivalent  integration 
scheme  becomes: 

T  /  /c ell  UdmM  +  Leu  ( FdO  -  Gdm )  =  [  icM  H^mdO  (3.1) 

Clt  J  Jarea  J faces  J  J  area 

where  U  is  the  state  vector,  F  and  G  are  the  convective  flux  vectors,  and  Hi  is  the 
inviscid  portion  of  the  source  term  vector  H.  Physically,  the  first  term  represents 
the  time  rate-of-change  of  the  state  vector  U  over  a  cell  with  area  A,  discretized  as 
The  subscript  c  indicates  the  cell  center  and  A Uc  =  —  U™  is  the  change 
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nw 


Figure  3.1:  Quadrilateral  cell  node  and  face  designations. 

in  Uc  over  a  single  time  step.  The  second  term  represents  the  contribution  of  the 
convective  fluxes  and  is  calculated  using  trapezoidal  integration.  For  example,  the 
convective  flux  on  the  west  face  of  the  cell  in  Fig.  3.1  is  given  by: 

+  J’*“)  (0„„  -  «.„)  -  (Gn"  +  —)  (m„w  -  m,w) 

where  the  subscripts  on  F  and  G  indicate  the  value  of  the  fluxes  at  these  corners. 
Finally,  the  third  term  represents  the  contribution  of  the  inviscid  source  term  over 
the  cell  area,  discretized  as  A HicA.  The  complete  discretization  of  Eqn.  3.1  is: 

+  {enw  -  e„)  -  ( A1+A1)  (m„„  _  m,ul) 

+  +2  f '-)  («„  - «„„) -  ) (m„ - m„) 

+  (f‘°  \  —)  («„  -  M  -  (g“  +  G“)  (m»»  -  m„t)  (3.2) 

+  (0,„  -  «„)  -  ( AdTfi)  (m„,  -  mse)  =  AHicA, 


26 


where  the  cell  area  A  is  given  as: 


A  —  ^  [(^rae  sui')  nw  9 se )  ij^nw  T7lsc}(@ne  ^sk;)]  •  (3.3) 

It  should  be  noted  that  A  as  given  above  is  not  a  true  area,  but  with  the  inclusion 
of  the  r  and  h  terms  in  the  state  and  flux  vectors,  forms  an  equivalent  cell  volume 
in  the  quasi-3D  coordinate  system. 

3.2  Viscous  Terms 

Considering  the  contributions  of  the  inviscid  and  viscous  portions  of  the  equations 
separately,  the  complete  change  in  the  dependent  variables  may  be  written: 

T,  //«„  Uimde  =  j.  f  L  U.imiS  +  J  I L  u,dmi),  (3.4) 

CLZ  J  Jarea  CLZ  J  Jarea  CLZ  J  Jarea 

where  the  two  terms  on  the  right  hand  side  of  Eqn.  3.4  represent  the  inviscid  and 
viscous  contributions  to  the  time  rate-of-change  of  the  state  vector  £/,  respectively. 
Again  using  a  stationary  grid  of  quadrilateral  cells,  the  integration  of  the  viscous 
terms  takes  the  form: 

T  /  [ecu  UvdmdO  =  Lell  ( Rd9  -  Sdm)  +  f  ftl  Hvdmd9.  (3.5) 

CLZ  J  Jarea  J  faces  J  Jarea 

Here  the  second  term  represents  the  contribution  of  the  diffusive  fluxes,  and  Hv  is 
the  viscous  portion  of  the  source  term. 

The  control  volume  used  for  the  viscous  terms  is  modified  slightly  from  that 
used  for  the  inviscid  terms.  The  aim  is  to  write  the  second  derivatives  in  the  viscous 
terms  in  such  a  manner  that  the  complete  second  derivative  at  each  node  point  can 
be  assembled  solely  from  the  information  available  in  the  surrounding  cells.  As  is 
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Figure  3.2:  Primary  cell  for  integration  viscous  terms. 

the  case  for  the  inviscid  terms  described  above,  these  derivatives  are  not  complete 
until  a  sweep  through  all  cells  has  been  accomplished. 

Considering  a  control  volume  around  point  o.  as  shown  in  Fig.  3.2,  the  line 
integration  about  the  cell  abed  may  be  written  (omitting  the  source  term  for  clarity): 
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=  {Rda^Oda  -  SdaAmda  +  RcdAOcd  -  ScdAmcd 

Aabcd 

“I-  Rbc^^bc  Sbc^'Tft'bc  ”F  -^ab^-^ab  Fq.6  A?7?,a5  J*  ?  (3.6) 

where  Rda  =  ~(Rd+Ra),  Sda  =  1(^+5^),  etc.  and  A0da  =  6^-0a,  A mda  =  md~ma , 
etc. 

In  order  to  rewrite  this  equation  in  terms  of  cell  variables,  the  following  geo¬ 
metric  terms  are  defined.  Considering  the  generic  cell  shown  in  Fig.  3.1,  with  the 
local  cell  grid  coordinates  j  and  k  as  shown,  define: 

Am j  —  ~{irfine  “h  mse  mnw  r^sw )  ? 

A 9j  —  ~(^ne  "F  $se  $niL>  $ sw)  •> 

Arrik  —  (^ne  “f~  ^nw  ^Y^se  sw )  9 

A6F  =  “(0ne  +  0nw  —  0se  —  y .  (3-7) 

Noting  that  A Qda  =  |(A0ii?+A0jc),  A0crf  =  |(A0fcc  +  A0fci?),  etc.,  Eqn.  3.6  becomes: 

At/„„  =  {  /UA0„  +  A0M)  -  Sia(/\m1D  +  Am,J 

+  Rcd(A6kc  +  A9kD )  —  Scd(  A/ kc  + 

—  Rbc(A9jg  +  A  9jc)  +  Sbc{ArrijB  +  A  rrijc) 

-  Rab(A9kA  +  A 9ks)  +  Sab(AmkA  +  A mkB)  }  .  (3.8) 

In  order  to  evaluate  the  values  of  R  and  S,  secondary  control  volumes  are 
drawn  about  the  midpoints  of  the  edges  of  the  primary  control  volume,  as  shown  in 
Fig.  3.3.  Considering  the  “east”  secondary  control  volume  in  particular  (Fig.  3.4), 
the  line  integration  for  a  first  derivative  term  becomes  (using  §|  as  an  example, 
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WEST 


SOUTH 


Figure  3.3:  Secondary  cells  for  viscous  term  integration, 
where  u  is  a  representative  velocity  component): 


du 


86 


cd 


1 

^ ad  dJ  cfbc 
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ad  d!  c1  be 


( u  dm) 

d'c'bc 


{ud(md'  —  mad )  +  u0'{mc'  -  md>) 


+uc(mbc  -  mci)  +  u0(mad  -  mbc )  }  . 


(3.9) 


If  the  assumption  is  made  that  Ac  ~  Ad  ~  Aadd’C'bc ,  then  with  the  definitions 
given  in  Eqn.  3.7,  Eqn.  3.9  may  be  split  into  two  pieces  that  only  contain  operations 


The  east  secondary  cell  contribution  to  AUVo  in  Eqn.  3.8  may  then  be  written  as: 


a  a 


At 


abed 


vo  I  east 


2^4  abed. 

Ate 


Rcd(  A6kc  +  A$kD) 
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2  A2C 
A to 
2AI 


\{u0’(mc’  -  m0i)  -  ucAm.jc  +  u0{m0  -  mbc)}  ZA6>fcc.] 
[{u0>(mo:  -  md>)  +  udAmjD  +  u0(mad  -  m0)}  A 0kD)  ■  (3.11) 


Similar  expressions  may  be  derived  for  the  S  term  in  the  east  secondary  cell  and  for 
all  terms  in  the  remaining  secondary  cells. 

Using  the  method  which  produced  Eqn.  3.11,  the  complete  expression  for  the 
viscous  fluxes  at  point  o  may  be  written  as  a  summation  of  four  terms,  one  for 
each  of  the  four  surrounding  cells.  Each  of  these  four  terms  will  use  information 
locally  available  to  that  cell  only;  no  information  from  neighboring  cells  is  required. 
Therefore,  the  viscous  terms  may  be  calculated  by  a  single  sweep  through  the  cells, 
with  the  contributions  from  each  secondary  cell  fragment  calculated  and  distributed 
to  the  appropriate  corner  nodes. 


3.3  Boundary  Conditions 

Boundary  conditions  for  the  solver  must  be  applied  at  the  inflow  and  outflow  of  the 
flow  domain,  along  solid  surfaces,  at  the  periodic  boundary  between  adjacent  blade 
passages,  and  at  the  sliding  interface  between  adjacent  blade  rows. 

3.3.1  Inflow  and  Outflow  Conditions 

The  boundary  conditions  specified  at  the  inlet  and  exit  of  the  domain  are  deter¬ 
mined  using  a  modification  of  the  characteristic  variable  approach  described  by 
Giles  [27,  28].  This  method  uses  a  knowledge  of  characteristic  theory  to  develop  the 
appropriate  boundary  conditions  for  the  flow  regime  under  consideration.  The  use  of 
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characteristic  theory  allows  the  development  of  non-reflective  boundary  conditions. 
These  are  important  in  turbomachinery  computations  due  to  the  close  proximity  of 
the  inlet  and  exit  boundaries  to  the  blade  surfaces. 

Two  different  sets  of  boundary  conditions  are  generally  used  in  turbomachinery 
computations.  Transonic  compressor  calculations  require  that  the  inlet  relative  Mach 
number  be  supersonic,  while  the  axial  Mach  number  remains  subsonic;  turbine  cal¬ 
culations  use  a  subsonic  inlet  Mach  number  condition.  Both  compressor  and  turbine 
cascade  computations  require  a  subsonic  exit  condition  with  a  specification  of  the 
downstream  static  pressure. 

A  complete  description  of  the  boundary  condition  development  and  implemen¬ 
tation  may  be  found  in  Appendix  B. 

3.3.2  Solid  Surface  Conditions 

On  solid  surfaces,  the  application  of  characteristic  theory  is  inappropriate  because 
the  assumption  of  inviscid  flow  is  no  longer  valid.  Viscous  calculations  require  that 
the  surface  velocity  of  the  fluid  match  that  of  the  wall  boundary,  the  so  called 
“no-slip”  condition.  While  this  condition  serves  to  specify  the  surface  velocities, 
the  surface  temperature  and  pressure  must  come  from  a  knowledge  of  the  phys¬ 
ical  boundary  conditions,  coupled  with  a  solution  of  the  momentum  equation  in  a 
direction  normal  to  the  surface.  Again,  compressor  and  turbine  cases  are  different. 
Compressor  airfoils,  which  are  generally  thin  and  operate  without  internal  cooling, 
are  assumed  to  use  an  adiabatic  surface  condition.  Turbine  airfoils,  on  the  other 
hand,  often  contain  internal  cooling  air  and  are  considered  to  operate  with  a  near 
constant  surface  temperature.  (A  full  thermal  analysis  of  the  coupled  fluid-structure 
domain  would  produce  an  improved  boundary  condition,  but  is  beyond  the  scope  of 
the  current  work). 
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With  the  surface  velocities  and  temperature  known,  a  solution  of  the  normal 
momentum  equation  provides  the  surface  pressure.  In  the  solution  of  this  equation, 
there  are  a  number  of  components  worth  discussing.  The  first  is  a  viscous  term 
that  remains  after  assuming  the  no-slip  condition  and  neglecting  streamwise  velocity 
gradients: 

dp_  _  \i_cN_  /3_12> 

drj  Re  dr}2 ' 

where  rj  is  the  direction  normal  to  the  wall.  The  derivative  of  the  velocity  component 
normal  to  the  wall,  u,  will  be  small  and  is  further  multiplied  by  the  inverse  of  the 
Reynolds  number.  Therefore,  this  term  is  typically  neglected  for  high  Reynolds 
number  flows. 

The  second  term  to  consider  is  the  pressure  gradient  resulting  from  wall  cur¬ 
vature.  In  this  case: 


dp 

drj 


r sy 


(3.13) 


where  u  is  the  streamwise  velocity  component  parallel  to  the  wall  and  Rc  is  the  local 
wall  radius  of  curvature.  This  term  might  become  important  for  the  small  radii  of 
typical  compressor  blade  leading  and  trailing  edges,  but  the  streamwise  grid  density 
required  to  resolve  the  boundary  layer  in  these  areas  will  usually  insure  that  this 
term  is  also  quite  small,  and  is  therefore  neglected  in  the  current  analysis. 

The  final  term  results  from  the  movement  of  the  wall,  but  only  for  cases  where 
the  radius  of  the  streamsheet  is  also  changing.  In  this  case  the  normal  momentum 
equation  yields: 

JT  =  -Pve  sin(aw)— ,  (3.14) 

drj  r 


where  ctw  is  the  wall  angle  measured  from  the  streamwise  direction.  For  a  stationary 
wall,  ve  =  0,  and  the  usual  condition  of  dp /dr}  =  0  is  obtained.  However  for  a  moving 
wall,  vg  =  fir,  and  this  term  can  no  longer  be  ignored.  Experience  has  shown  that, 
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Figure  3.5:  Treatment  of  periodicity  condition. 

for  a  typical  compressor  rotor  computation,  this  term  results  in  a  1-2%  increase  in 
the  surface  pressure  from  the  dp/ dp  —  0  condition. 

3.3.3  Periodicity 

In  order  to  simulate  an  infinitely  repeating  cascade,  a  periodic  boundary  condition 
is  implemented  upstream  of  the  leading  edge  and  downstream  of  the  trailing  edge 
of  the  blade.  Nodes  that  lie  along  these  boundaries  are  treated  as  internal  points  in 
the  integration,  but  an  extra  pass  is  required  to  combine  fluxes  across  the  periodic 
boundary  and  equate  the  values  at  both  points.  Considering  Fig.  3.5,  the  fluxes 
calculated  in  cells  A  and  B  during  the  normal  flux  calculation  pass  are  sent  to  node 
o,  those  in  cells  C  and  D  are  sent  to  d .  The  second  pass  adds  the  partial  fluxes  at 
nodes  o  and  d  and  sets  both  nodal  values  to  this  correct  total  flux.  This  double  pass 
is  performed  during  both  inviscid  and  viscous  flux  evaluations,  as  well  as  during  the 
calculation  of  the  artificial  dissipation  terms. 
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For  cases  which  contain  more  than  one  passage  in  a  blade  row,  the  periodicity 
condition  is  still  applied,  but  only  along  the  top  and  bottom  boundaries  of  the 
multiple  passage  grid  block.  The  grid  lines  upstream  and  downstream  of  the  blades 
in  the  middle  of  the  multiple  passage  block  are  implemented  as  internal  points,  and 
no  special  treatment  is  required. 

3.3.4  Intra-Blade  Row  Interface  Conditions 

The  interface  between  adjacent  blade  rows  occurs  along  an  overlap  region  of  two 
grid  lines,  where  the  overlapping  lines  lie  directly  on  top  of  one  another.  Requiring 
the  grid  lines  to  slide  on  top  of  one  another  significantly  simplifies  the  interpolation 
procedure,  which  may  be  implemented  as  a  series  of  line  interpolations  rather  than 
a  more  expensive  area  interpolation.  The  speed  at  which  the  interpolation  is  done 
is  offset  by  the  fact  that  it  must  be  done  twice  in  each  overlap  region  of  a  blade 
row  for  each  stage  of  the  n-stage  Runge-Kutta  time  integration  procedure.  Since  all 
blade  rows  in  the  calculation  are  computed  in  the  absolute  reference  frame,  there  is 
no  need  to  impose  velocity  triangle  relationships  at  the  row  interfaces. 

The  treatment  of  the  flux  variables  at  the  interface  is  done  in  a  conservative 
manner  similar  to  that  of  Rai  [52,  53].  The  procedure  simply  requires  an  integration 
of  the  flow  variables  before  and  after  the  interpolation  between  grids.  The  inter¬ 
polated  values  are  then  adjusted  by  the  ratio  of  the  two  integrated  values.  This 
insures  that  both  the  distribution  of  the  variable,  as  well  as  its  integrated  value,  are 
preserved  across  the  interface  boundary. 

Choice  of  the  interpolation  algorithm  is  important,  particularly  when  the  grid 
is  relatively  coarse.  Not  only  is  it  important  to  conserve  the  flow  variables  across 
the  interface,  but  it  is  also  necessary  to  pass  their  profiles  with  as  little  distortion 
as  possible.  To  this  end,  three  different  interpolation  functions  were  implemented, 
each  based  on  Lagrangian  polynomials,  differing  only  in  the  number  of  points  used 
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in  the  interpolation.  In  general,  the  value  of  a  function  at  some  point  f(x)  may 
be  found  from  values  at  surrounding  points,  /(x0),  f(xx), . . .  by  using  the 

following  formula  [14,  pg.  27]: 


N 

f(x )  =  E 


i= 0 


N 

n 

j— ° 


0  z 

( Xi  Xj ) 


N 


f(xi )  =  I 2Li(x)f(xi), 
1=0 


(3.15) 


where  Li(x)  are  the  Lagrange  polynomials.  For  example,  a  three  point  interpolation 
is  written  as: 


(^0  Zi)(z0  ^2)  (^1  ^2)(^1  ^0) 

(x  -  x0)(x  -  xx) 


(x2  -  X0){x2  -  xx) 


f(x 2)- 


(3.16) 


This  method  of  interpolation  uses  a  polynomial  of  degree  N  to  locally  fit  the  existing 
data  and  determine  the  value  of  the  desired  point;  i.e. ,  N  =  1  uses  a  linear  fit,  N  =  2 
uses  a  quadratic  function,  etc.  Interpolation  has  been  implemented  with  N  =  1,2, 
and  3  corresponding  to  two,  three,  and  four  points  respectively.  In  each  case,  the 
points  used  in  the  interpolation  are  found  by  requiring  that  the  desired  point  lie 
between  the  two  center  points;  the  third  point  in  a  three-point  interpolation  is  located 
at  the  closest  of  the  next  two  outer  points.  The  determination  of  the  Lagrange 
polynomials  is  made  once  each  time  the  grid  position  is  indexed;  these  values  are 
stored  and  reused  for  each  interpolation  performed  at  this  grid  position. 

An  example  of  the  value  of  the  higher-order  interpolation  may  be  found  in  the 
density  profile  from  the  hot-streak  calculations  described  in  Sect.  6.1.4.  Starting 
with  identical  profiles,  an  interpolation  is  performed  using  each  of  the  Lagrangian 
schemes;  the  resulting  profiles  are  shown  in  Fig.  3.6.  Not  surprisingly,  the  three-point 
algorithm  provides  a  significantly  better  profile  near  the  center  of  the  profile  while 


all  methods  are  indistinguishable  in  the  linear  regions.  The  four-  and  three-point 
interpolations  are  nearly  identical,  as  expected,  because  of  the  quadratic  nature  of 
the  profile.  Since  the  flow  equations  are  nonlinear,  and  are  at  most  of  quadratic 
order,  the  three-point  function  should  be  sufficient  to  predict  most  smooth  profiles, 
and  in  practice  requires  little  computational  overhead  beyond  that  of  the  two-point 
function. 


3.4  Artificial  Dissipation 

The  explicit  addition  of  artificial  dissipation  is  generally  required  for  time-marching 
schemes  with  spatial  central  differencing.  Damping  is  needed  both  to  suppress  spu¬ 
rious  oscillations  in  smooth  regions  of  the  flow,  as  well  as  to  capture  flow  disconti¬ 
nuities,  such  as  shocks.  This  section  examines  the  method  in  which  the  second-  and 
fourth-order  damping  terms  are  added  to  the  code. 

Second-order  smoothing  is  required  by  central-difference  schemes  in  order  to 
capture  discontinuities  in  the  flow,  such  as  shock  waves.  While  a  shock  would  ideally 
have  zero  thickness,  central  difference  schemes  do  not  permit  this  type  of  discon¬ 
tinuity  as  a  valid  solution.  The  introduction  of  second-order  smoothing  serves  to 
locally  change  the  difference  equations  to  first-order,  which  do  permit  discontinuities. 
Shocks  are  smeared  over  several  grid  points  with  a  minimum  amount  of  oscillation 
(undershoot  and  overshoot)  up-  and  downstream  of  the  shock. 

Fourth-order  smoothing  is  applied  to  suppress  the  formation  of  the  odd-even,  or 
“sawtooth,”  modes  characteristic  of  central  difference  schemes.  A  central-difference 
scheme  uses  a  three-point  stencil  in  order  to  evaluate  first  derivatives,  but  does  not 
use  the  value  at  the  center  node.  This  causes  an  uncoupling  of  values  at  alternate  grid 
points,  resulting  in  the  formation  of  the  so  called  odd-even  modes  shown  in  Fig.  3.7. 
Using  a  three-point  stencil  to  calculate  the  first  derivative  “■  of  the  solution  in  this 
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Figure  3.6:  Comparison  of  interpolation  polynomials  in  hot  streak  calculation 


Figure  3.7:  Odd-even  or  sawtooth  mode. 

figure  would  yield  a  value  of  zero  at  all  points.  Therefore,  this  solution  would  be 
accepted  as  a  valid  numerical  solution  of  |£  =  0,  while  being  wholly  unsatisfactory. 
Fourth-order  damping  uses  a  five-point  computational  stencil  to  detect  and  suppress 
these  oscillations. 

The  smoothing  used  in  the  current  effort  is  an  unstructured  implementation 
of  the  damping  terms  normally  found  in  finite  difference  solutions.  The  smoothing 
terms  are  explicitly  added  to  the  right  hand  side  of  Eqn.  2.15.  The  combined  second- 
and  fourth-order  smoothing  operator  has  the  form: 

D(U)  =  Dj(U)  +  Dk(U),  (3.17) 

where  again  j  and  k  are  the  cell  local  grid  directions.  In  the  ^-direction,  the  operator 
is  given  by: 

Dj  ( U )  -  ft  (e26]U  -  e4S*u)  ,  (3.18) 

where  the  6?  and  Sj  operators  represent  the  second  and  fourth  differences  in  the 
ji-direction.  A  similar  expression  may  be  written  for  the  ^-direction.  In  practice,  a 
sweep  is  made  through  the  computational  domain  to  calculate  the  second  differences 
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at  all  nodes.  The  fourth  difference  is  then  found  with  another  sweep  to  “second 
difference  the  second  differences.”  The  coefficients,  /3j  and  /?*.  are  approximations 
of  the  cell  spectral  radius,  scaling  the  added  dissipation  to  the  change  in  the  flux 
residual.  One  approach  is  use  the  values  to  cancel  the  At /A  term  which  results 
from  the  integration  of  Eqn.  2.15.  In  this  case,  the  /3’s  have  the  form  [36,  39]: 

ft  =  A  =  -A-  (3-19) 


Kroll  et  al.  [39]  define  At*  as  the  local  cell  time  step  based  on  an  unit  Courant- 
Freidrichs-Lewy  (CFL)  number,  where  the  CFL  number  is  defined  as  the  ratio  of 
the  time  step  of  the  numerical  scheme  to  the  transit  time  across  a  characteristic 
cell  dimension  at  the  local  maximum  wavespeed  [17].  An  equivalent  form  is  given 
by  Jorgenson  and  Chima  [36],  who  use  the  local  time  step  calculated  as  a  function 
of  the  stability  limit  of  the  cell,  but  use  a  subsequent  increase  in  the  second-  and 
fourth-order  damping  coefficients  described  below.  An  alternate  form  is  given  by 
Mavriplis  and  Jameson  [47],  who  scale  the  dissipation  depending  on  the  stretching 
in  the  local  grid  coordinates.  A  modified  form  of  this  latter  approach  is  used  in  the 
current  effort,  with: 


Pi  —  + 

Pk  —  Afc  ( 1  + 


(3.20) 


Here,  A j  and  A*,  are  the  maximum  eigenvalues  in  the  two  local  grid  directions: 
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where  c  is  the  local  speed  of  sound,  and  A 0j,  Arrij,  etc.  are  given  by  Eqn.  3.7. 

The  62  and  64  coefficients  are  given  by: 

e2  =  c2AP,  (3.22) 

e4  =  max(0,  a4  —  e2).  (3.23) 

The  pressure  switch  A P  is  implemented  as  a  normalized  second  difference  in  pressure 
at  a  particular  node: 

Pn') 

A  Pi  = - 

E*#n(P»-  +Pn) 

where  the  n  points  are  the  edge  nodes  around  the  four  cells  bordering  a  particular 
node.  Near  a  shock  wave,  where  the  fourth-order  smoothing  becomes  destabilizing, 
the  pressure  gradient,  and  therefore  A P,  becomes  large.  This  causes  the  <t2AP  term 
to  dominate  <74  in  Eqn.  3.23,  switching  off  the  fourth-order  smoothing.  In  smooth 
regions  of  the  flow,  A P  is  small.  The  dissipative  terms  are  of  third  order  and  do 
not  reduce  the  accuracy  of  the  central-difference  scheme.  Near  shock  waves,  the 
pressure  switch  becomes  large  and  the  dissipation  becomes  locally  of  first  order. 
The  coefficients  cr2  and  a4  are: 

n  =  0(1). 

cr4  =  0(i).  (3.25) 

The  smoothing  terms  are  further  modified  by  using  a  velocity  scaling  in  order 
to  avoid  contamination  of  the  viscous  terms  in  the  boundary  layer.  Both  smoothing 
terms  are  multiplied  by  a  factor  which  scales  with  the  ratio  of  the  local  velocity 
to  that  determined  from  the  local  static-to-total  pressure  ratio  [21].  For  rotational 
flows,  care  must  be  taken  to  insure  that  this  ratio  is  formed  with  the  local  relative 


(3.24) 
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velocity  and  total  pressure  to  insure  proper  detection  of  the  boundary  layer  on  a 
rotating  blade. 


3.5  Temporal  Integration  Scheme 

The  following  sections  describe  the  numerical  integration  scheme  used  to  advance  the 
flow  equations  forward  in  time.  The  first  section  describes  the  dual  time  stepping 
algorithm,  which  allows  an  explicit  multi-stage  time-stepping  method  to  be  used 
as  a  driver  for  a  fully  implicit  scheme.  A  description  is  then  presented  of  the  two 
convergence  acceleration  techniques  implemented  in  the  current  effort. 

3.5.1  Dual  Time  Stepping 

While  explicit  algorithms,  and  particularly  multi-stage  Runge-Kutta  schemes,  have 
a  low  operation  count  per  iteration,  they  suffer  a  serious  time  step  restriction  as 
compared  to  their  implicit  counterparts.  This  condition  is  further  exacerbated  in 
unsteady  flows  where  the  global  minimum  time  step  must  be  used  throughout  the 
domain.  In  these  cases,  there  may  be  a  variation  of  several  orders-of- magnitude  in  the 
time  step  required  for  stability  between  the  largest  and  smallest  cells.  Additionally,  a 
direct  temporal  integration  of  Eqn.  2.15  prohibits  the  use  of  convergence  acceleration 
techniques,  such  as  local  time  stepping,  implicit  residual  smoothing,  or  multigrid, 
for  unsteady  flows. 

To  overcome  these  limitations,  the  Navier-Stokes  equations  may  be  reformu¬ 
lated  using  one  of  the  dual  time  step  methods  described  by  Jameson  [32]  and  Weiss 
and  Smith  [76].  These  methods  use  a  standard  explicit  scheme  as  the  driver  for  a 
fully  implicit  time  stepping  scheme,  differing  only  in  the  details  of  their  implementa¬ 
tion  (which  will  be  discussed  in  Section  3.5.3  below).  Updates  of  the  solution  domain 
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are  performed  at  time  intervals  determined  by  the  required  temporal  accuracy  of  the 
solution,  rather  than  those  prescribed  by  the  smallest  cells  in  the  domain. 

A  compact,  differential  form  of  Eqn.  2.15  may  be  written  as: 

A^  +  TZ(U)  =  0,  (3.26) 

where  7 Z(U)  is  the  residual  containing  both  the  convective  and  diffusive  fluxes,  as 
well  as  the  artificial  dissipation  terms.  In  a  steady  flow  case,  both  terms  in  a  dis¬ 
cretized  form  of  Eqn.  3.26  tend  toward  zero  as  the  solution  converges.  Acceleration 
methods  are  effective  because  they  drive  these  terms  to  zero  with  fewer  iterations, 
leaving  no  residual  effect  on  the  solution.  However,  transients  computed  using  these 
acceleration  methods  are  completely  non-physical  since  time  has  not  been  advanced 
uniformly  across  the  mesh.  Therefore,  when  using  these  techniques,  the  variable  t 
in  Eqn.  3.26  must  be  considered  solely  as  an  integration  variable  and  not  as  physical 
time.  This  point  will  become  important  in  the  construction  of  the  dual  time  step 
formulation. 

In  an  unsteady  flow  case,  neither  term  in  Eqn.  3.26  is,  in  general,  zero,  yet 
the  sum  of  the  two  terms  is  still  zero  for  all  time.  For  temporal  accuracy,  time 
must  be  advanced  uniformly  throughout  the  mesh,  disallowing  the  use  of  any  form 
of  acceleration  technique.  In  the  construction  of  a  dual  time  step  method,  a  new 
residual  is  defined  which  includes  both  terms  in  Eqn.  3.26.  This  combined  residual 
is  then  driven  to  zero  using  a  multi-stage  Runge-Kutta  algorithm.  Recasting  the 
equations  in  this  way  produces  a  formulation  analogous  to  that  used  in  the  solution 
of  the  steady-state  equations,  and  again  allows  the  use  of  acceleration  methods. 

A  fully  implicit  formulation  of  Eqn.  3.26  may  be  approximated  as: 

ADt{U{n+1) )  +  H(U{n+1))  =  0,  (3.27) 
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where  the  superscript  denotes  the  time  level  and  Dt(U )  is  a  backwards-difference 
temporal  operator  of  some  specified  accuracy.  Using  a  second-order  backwards- 
difference  operator,  Eqn.  3.27  becomes: 


A 

2A  t 


3£/("+l)  _  4 jj{n)  _|_  jj(n- 1) 


+  R(t/(*+1>)  =  0, 


(3.28) 


which  is  the  combined  residual  previously  described.  Because  this  combined  residual 
will  be  driven  to  zero  in  a  converged  solution,  Eqn.  3.28  may  be  treated  as  a  modified 
steady-state  problem,  solved  by  stepping  in  a  fictitious  time  r: 


-  A 


du(n+v 

d  r 


A 

2A  t 


3 U(n+1)  -  4U(n)  +  U{n~l)]  +  K(U{n+1))  =  7l*(t/(n+1))  (3.29) 


or 

flTj(n+l) 

A — - - +  7T(U(”+1))  =  0,  (3.30) 

Or 

which  is  of  the  same  form  as  Eqn.  3.26.  Between  each  increment  in  real  time  the 
solution  converges  in  the  fictitious  time  variable  r,  causing  the  dU^71^/ dr  term  to 
vanish,  resulting  in  a  proper  solution  of  Eqn.  3.28.  Selection  of  the  physical  time 
increment  At  is  now  driven  by  the  desired  time  accuracy  of  the  solution,  rather  than 
by  the  global  minimum  time  step  restriction.  The  progression  of  the  solution  in  ficti¬ 
tious  time  does  not  affect  the  solution  in  real  time,  therefore  any  of  the  conventional 
convergence  acceleration  techniques  may  be  used.  In  the  current  implementation, 
both  local  time  stepping  and  implicit  residual  smoothing  (described  in  the  following 
section)  have  been  implemented. 

Incorporating  the  dual  time  stepping  into  a  standard  Runge-Kutta  code 
requires  little  additional  effort.  The  backward-difference  operator  of  Eqn.  3.29  is 
included  as  an  additional  source  term.  This  source  term  and  the  storage  of  the  solu¬ 
tion  at  the  previous  time  level,  are  the  only  modifications  required.  The  complete 
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algorithm  for  the  dual-time  iterations  is  outlined  in  Fig.  3.8,  where  the  number 
of  iterations  in  real-time  is  specified  with  the  nreal  parameter  and  nfake  is  the 
number  of  pseudo-time  iterations  required  to  converge  to  the  solution  of  U^n+l\ 

Initialization  of  the  value  of  U^n+1')  at  the  beginning  of  a  sequence  of  ficti¬ 
tious  time  steps  is  performed  by  predicting  its  value  using  the  previous  value  of  the 
time  derivative.  In  practice,  it  has  been  found  necessary  to  limit  the  change  in  the 
value  of  a  dependent  variable  to  approximately  30%  of  its  current  value  to  prevent 
destabilizing  oscillations  of  flow  variables  which  have  values  near  zero. 

Use  of  the  dual  time  step  integration  scheme  requires  the  specification  of  a 
suitable  real  time  step  which  is  directly  related  to  the  desired  temporal  accuracy 
of  the  solution.  For  periodic  solutions,  the  time  step  is  specified  by  a  division  of 
the  period  of  the  solution  into  a  number  of  equal  increments;  non-periodic  solutions 
require  a  direct  specification  of  the  real  time  step  value.  The  choice  of  the  real  time 
step  is  a  topic  to  be  discussed  in  Chapter  6. 


3.5.2  Convergence  Acceleration 

Two  techniques  for  accelerating  the  convergence  of  the  numerical  scheme:  local  time 
stepping  and  implicit  residual  smoothing.  Both  of  these  methods  are  applicable 
to  flow  cases  converging  to  a  steady-state  or  those  using  the  dual  time  stepping 
described  in  the  previous  section. 

Local  time  stepping  allows  each  cell  to  use  its  maximum  permissible  time  step 
based  on  a  local  stability  analysis.  Viscous  flow  cases  require  that  both  the  convective 
and  diffusive  terms  be  included  in  this  analysis.  Assuming  the  use  of  the  dual  time 
stepping  algorithm,  the  local  pseudo-time  step  is  given  by: 


At  =  CFL 


/  A  Tj  A  Tv  \ 
V  ATj  +  A Ty  / 


(3.31) 
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initialize 
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Figure  3.8:  Dual  time  stepping  algorithm. 


where  CFL  is  the  allowable  Courant-Friedrichs-Lewy  stability  limit  of  the  particular 
time  integration  scheme,  and  At;  and  A tv  are  the  local  convective  and  diffusive  time 
step  limits,  respectively.  The  convective  limit  used  here  is  similar  to  that  described 
by  Jorgenson  and  Chima  [36]  for  calculation  of  quasi-3D  flows,  with  modifications 
into  the  appropriate  finite  volume  form: 


At; 


A 


A;T  I  wo  I  A k  +  cF A?  +  A?  + 


vgdr  /  dm  A2 
2c, /A?  +  A  l 


where 


A,  = 

Afc  = 


I  Am,- 

r 

1  A 

r 


I  Mj 


a  ek 


(3.32) 


and  c  is  the  local  sound  speed.  The  viscous  limit  is  a  finite  volume  formulation  of 
the  diffusive  limit  given  by  Arnone,  et  al.[4]: 


A  Tv  — 


A 2 

2.5y  /  n  /ir  \ 

p  \Pr  Prx ) 


(3.33) 


All  quantities  are  assumed  to  be  cell-averaged  values. 

Implicit  residual  smoothing  (IRS)  is  used  to  remove  spikes  from  the  flux  resid¬ 
uals,  reducing  the  high-frequency  errors.  Use  of  IRS  can  expand  the  stability  region 
and  improve  the  damping  characteristics  of  multi-stage  time  stepping  procedures. 
For  two-dimensional  problems,  the  total  flux  residual  7 Z*(U)  of  Eqn.  3.29  is  replaced 
by  a  solution  of  the  implicit  equations  [47]: 


n*  =  1l*  +  eV2^ 


(3.34) 
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where  i  is  the  node  point  index,  V277*  is  the  undivided  Laplacian  of  the  residual, 
and  e  is  the  smoothing  coefficient.  While  the  smoothing  can  be  applied  using  a 
uniform  value  of  e,  coefficients  which  are  functions  of  the  local  grid  stretching  and 
flow  conditions  have  been  shown  to  greatly  improve  the  rate  of  convergence  [51]. 
The  use  of  locally  varying  coefficients  makes  the  scheme  more  implicit  in  the  grid 
direction  of  higher  stretching,  and  less  implicit  in  the  direction  with  less  stretching. 
In  an  unstructured  quadrilateral  mesh,  the  Laplacian  operator  of  Eqn.  3.34  may  be 
replaced  by: 

+  ej  8fR*  +  8fR~,  (3.35) 

where  again  j  and  k  are  the  characteristic  directions  within  the  cell  and  82  is  the 
second  difference  operator  in  the  appropriate  direction.  The  smoothing  coefficients 
are  [47]: 


1  f  (  CFL 


&  - 1  ,  0  , 


62  =  maX  [4  [CFL*P> 


(  1  /  CFL  0  V  ,  \ 

ek  =  max  -  77 7777  Pk  )  ~  1  ,  0  , 


4  \\CFL* 


where  and  3k  are  the  scaling  coefficients  given  by  Eqn.  3.20.  The  7777  term  is 
the  ratio  of  the  CFL  numbers  in  the  smoothed  and  unsmoothed  schemes  and  has 
been  found  to  have  an  optimum  value  of  ~  2  [72].  The  second  difference 

operators  are  similar  to  those  used  in  the  artificial  dissipation  and  are  formed  in 
an  analogous  manner.  Two  Jacobi  iterations  are  used  to  provide  an  approximate 
solution  to  Eqn.  3.35;  the  smoothing  operator  is  applied  during  each  stage  of  the 
Runge-Kutta  procedure.  In  practice,  the  use  of  IRS  technique  reduces  the  number 
of  pseudo-time  iterations  required  to  converge  the  solution  at  each  real  time  step 
by  approximately  one-third,  with  only  a  small  increase  in  the  computational  time 
required  for  each  iteration. 
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3.5.3  Numerical  Scheme 


The  time  integration  algorithm  used  in  the  present  work  is  an  n-stage  Runge-Kutta 
method.  In  the  classical  Runge-Kutta  scheme,  the  residuals  are  evaluated  and  stored 
at  each  stage,  then  combined  in  the  last  stage  to  update  the  solution.  Storage  of 
these  residuals  can  be  quite  expensive  and  normally  leads  to  the  use  of  a  modified 
Runge-Kutta  scheme  where  a  given  stage’s  residuals  are  used  only  in  the  evaluation 
of  the  next  stage. 

The  Runge-Kutta  scheme  used  to  update  the  solution  in  pseudo-time  is  one  of 
the  modified  schemes  of  Jameson  [32].  In  general,  an  r-stage  scheme  may  be  given 
by: 

_  jj{n) 

Lr(n+l,l)  =  jj{n)  _ai^L  1q(0)  +  £(0)1 


U{n+ 1’9)  =  c/(n)  -  aq-j-  [Q(<?_1)  + 
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D{q)  =  pqD  [U{n+1’9)]  +(1  -/39)D(9-1).  (3.41) 

Here  Q  [17]  is  the  combined  inviscid  and  viscous  flux  residual  and  D  [U]  is  the  artificial 
dissipation  operator.  Two  modified  schemes  are  given  by  Jameson  [32].  The  first  is 
a  four-stage  scheme  with  two  evaluations  of  the  artificial  dissipation: 


1 

0l  =  3’ 

/? i  =  l 

4 

Oft  =  - . 

15' 

*  =  \ 

5 

ft  =  0 

=  1, 

O 

II 

(3.42) 

The  second  is  a  five-stage  scheme  with  three  dissipation  evaluations: 

^i  =  l 

«2  =  &  =  0 

6 

a3  =  /?3  =  0.56 

O 

«4  =  &  =  0 

a5  =  1,  &  =  0.44.  (3.43) 

Both  of  these  schemes  have  first  order  temporal  accuracy  (in  pseudo-time)  and 
yield  a  maximum  CFL  =  3.0  for  the  four-stage  scheme,  and  CFL  =  4.0  for  the 
five-stage  scheme.  These  values  are  typically  reduced  to  2.8  and  3.5  respectively 
with  the  inclusion  of  added  dissipation  [36,  39].  Both  of  these  schemes  have  been 
used  in  the  current  effort:  the  five-stage  is  commonly  used  for  the  flow  equations 
and  for  the  one-equation  turbulence  model  on  coarse  grids.  Adapted  meshes  often 
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require  the  use  of  sub-iterations  of  the  four-stage  scheme  for  the  integration  of  the 
one-equation  turbulence  model  (discussed  in  Chapter  4). 

In  Jameson’s  implementation  of  the  dual  time  stepping  method  [32],  the  Runge- 
Kutta  scheme  to  solve  Eqn.  3.30  in  the  form  as  described  above.  The  temporal 
derivative  source  term  is  evaluated  with  the  tTn+1  1  from  the  previous  stage.  The 
evaluation  of  the  q-th  stage  is  given  by: 


U(n+hg)  =  U{n)  -  aq— -  (Q{q~1]  +  D{q~l) 

/x  ' 

_  —  \oTT(n+l,q-l)  _  4 Tj(n)  ,  ly(n-l)]  \ 

2 At  l  \) 


Arnone,  et  al.[4]  have  noted  that,  although  this  time  discretization  of  Eqn.  3.29 
is  implicit,  stability  problems  can  occur  in  Jameson’s  method  if  the  time  step  in 
pseudo-time  exceeds  the  physical  time  step.  Therefore,  the  time  step  r  must  be 
limited  [4]: 

(3.45) 

for  solutions  without  multigrid.  The  effect  of  this  limitation  normally  only  appears 
in  viscous  calculations,  where  the  time  steps  in  the  inviscid  core  region  of  the  flow 
can  often  be  one  to  two  orders-of-magnitude  higher  than  those  in  the  viscous  layer. 

An  alternate  implementation  of  dual  time  stepping  is  given  by  Weiss  and  Smith 
[76]  who  uses  a  partially  implicit  treatment  of  the  temporal  derivative  term.  In  this 
case,  the  q-th  stage  of  the  Runge-Kutta  scheme  becomes: 


3Ar\ 
2A  t) 


A  U 


-aq~(Q^  +  D^ 


(3.46) 
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where  A U  =  U^n+1,q^  —  U^.  This  scheme  does  not  require  the  pseudo-time  step 
limitation  described  above. 

An  alternate  scheme  has  been  used  for  the  calculations  presented  in  this  thesis. 
This  method  is  a  simple  reformulation  of  Eqn.  3.44  with  a  fully  implicit  treatment 
of  the  temporal  derivative: 

(J+  3°AfT)  f/(”+‘:,)  =  f/(">  “  +  -D<,~1> 

+  At  K”>-Ct”-')]).  (3.4T) 

Like  the  scheme  of  Weiss  and  Smith,  this  scheme  has  not  exhibited  the  instabilities 
which  require  the  limitation  on  the  pseudo-time  step. 

While  all  of  these  schemes  have  been  implemented  in  the  computer  code  devel¬ 
oped  for  this  thesis,  studies  performed  to  determine  the  adequate  level  of  temporal 
resolution  produced  spurious  behavior  with  Jameson’s  method.  In  particular,  solu¬ 
tions  with  increasingly  smaller  real  time  steps  should  asymptotically  approach  the 
solution  obtained  using  a  standard  explicit  formulation.  In  some  cases,  however, 
refining  the  time  step  produced  an  increasing  level  of  divergence  of  the  dual  time 
step  solution.  This  phenomena  remains  unexplained;  it  has  never  been  observed 
when  using  either  the  formulation  of  Weiss  and  Smith  or  that  in  Eqn.  3.47. 

3.6  Properties  of  the  Unstructured  Scheme 

The  unstructured  time  integration  scheme  is  now  examined  with  respect  to  its  sta¬ 
bility,  accuracy,  and  conservation. 
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3.6.1  Stability 


In  order  to  examine  the  stability  properties  of  the  multi-stage  Runge-Kutta  time 
integration  scheme,  the  following  model  problem  is  considered: 

Ut  +  ux  jiAx  uxxxx  —  6,  (3.48) 

where  the  pAx3uxxxx  term  represents  an  added  third-order  dissipation.  With  fi  =  0, 
this  equation  represents  the  standard  one-dimensional  wave  equation  with  unit  char¬ 
acteristic  speed.  If  the  spatial  derivatives  are  approximated  with  central  differences, 
then  the  following  system  of  first-order  differential  equations  is  obtained: 

—Uj  +  Pj  =  0,  (3.49) 


where 


A*  -Pj  -  —  (uj+i-Wj-i) 

-  p^~(uJ+ 2  —  4uJ+i  +  6 uj  -  4 Uj-i  +  Uj- 2)  (3.50) 

and  A  =  At/ Ax  is  the  one-dimensional  CFL  number. 

Considering  the  von  Neumann  stability  analysis  presented  in  [39],  the  solution 
to  Eqn.  3.50  may  be  represented  by  the  summation  of  a  series  of  Fourier  modes 
expanded  in  the  x-direction: 


u  ^Y.Ukeipk{kAx\,  (3.51) 

k 

where  Uk  is  the  amplitude,  pk  is  the  wave  number  of  the  k-ih  mode,  and  the  number 
of  terms  k  is  equal  to  the  number  of  grid  points  in  the  mesh.  As  this  is  a  linear 
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equation,  the  effect  of  a  single  Fourier  mode  may  be  considered  and  the  final  solution 
obtained  through  linear  supposition.  Substituting  a  single  mode  into  Eqn.  3.49, 
yields: 


A tjt  (Ueip(jAx))  =  -At  •  P  (Ueip(jAx)) 

=  z{£)Ueip{jAx),  (3.52) 

where  £  =  pAx  is  the  phase  angle  and  z(£)  is  the  Fourier  symbol  of  the  spatial 
discretization  given  by: 


*(0 


-A 


4yu(l  —  cos  £)2  +  i  sin  £ 


(3.53) 


The  amplification  factor  of  the  time  integration  scheme  is  then  defined  as: 


(jn+l 
JJn  ’ 


(3.54) 


and  the  stability  region  is  defined  as  the  region  in  the  complex  plane  where  \g(z)\  < 
1.  For  the  scheme  to  be  stable,  z(£)  must  lie  within  this  stability  region  for  all 
— 7r  <  £  <  7r.  A  thorough  analysis  and  description  of  the  stability  boundaries  for 
various  Runge-Kutta  schemes  is  presented  in  [39]. 


3.6.2  Accuracy 

The  basic  central  difference  scheme  has  second-order  spatial  accuracy  on  uniform 
meshes  [37].  However,  grid  non- uniformities,  particularly  grid  stretching  and  grid 
skew,  can  lower  the  local  accuracy  to  first-order.  In  both  the  inviscid  and  viscous 
portions  of  the  flow  domain,  the  integration  of  flow  variables  to  form  either  first 
or  second  derivatives  is  performed  around  the  boundaries  of  some  computational 
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cell.  It  is  assumed  in  this  integration  process  that  the  center  of  the  cell  lies  at  the 
node  point  at  which  the  quantity  is  desired.  Green’s  function  is  used  to  evaluate 
these  derivatives,  yielding  a  cell-average  derivative  value  at  the  true  centroid  of  the 
cell.  Grid  stretching  reduces  the  accuracy  by  shifting  the  centroid  of  the  cell  away 
from  the  corresponding  node  point.  This  may  be  seen  by  performing  a  Taylor’s 
series  expansion  of  the  first  derivative  about  some  point  (x0,  y0).  Expanding  Green’s 
function: 


j  J  L^{x°’v°)dA 

1  f  f  d^u 

AJjJX-X°WX°’y°)clA 

U{x°’y°)dA+-'-  (3-55) 

For  all  of  the  higher  order  terms  to  cancel  we  must  have 

J  J  (x  -  x0)m(y  -  y0)ndA  =  0  (3.56) 

for  all  m,n  =  0,1,2,...,  implying  that  (x,y)  must  lie  at  (x0,y0).  In  other  words, 
the  use  of  Green’s  function,  when  the  true  cell  center  does  not  correspond  to  the 
assumed  nodal  point  center,  will  introduce  the  higher-order  error  terms  of  Eqn.  3.55. 
This  shift  in  the  node  point  from  the  true  cell  center  arises  through  stretching  in 
the  grid  as  shown  in  Fig.  3.9.  The  viscous  portion  of  the  solver  is  particularly 
sensitive  to  grid  stretching  errors  due  to  the  use  of  the  secondary  cells.  In  practice,  a 
maximum  stretching  factor  of  1.2  has  been  found  to  yield  acceptable  results  in  regions 
where  viscous  terms  are  important,  with  considerably  larger  stretching  allowed  in 
the  inviscid  portion  of  the  domain. 


du 
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Figure  3.9:  Shift  of  cell  centroid  through  grid  stretching. 

As  will  be  discussed  in  Chapter  5,  mesh  adaptation  can  have  an  additional 
impact  on  solution  accuracy,  but  careful  selection  of  the  divided  cell  geometry  can 
aid  in  an  overall  reduction  of  grid  stretching  errors. 

3.6.3  Conservation 

In  addition  to  the  consideration  of  a  scheme’s  accuracy,  consistency,  and  stability, 
one  must  also  insure  that  the  dependent  variables  are  globally  conserved  within 
the  computational  domain.  Demonstrating  this  will  show  a  scheme’s  conservation 
property .  The  Navier-Stokes  equations  express  conservation  of  mass,  momentum, 
and  energy.  To  show  global  conservation  in  a  discrete  solution  of  these  equations, 
the  equations  must  satisfy  an  equivalent  discrete  equation  that  approximates  the 
integral  equation.  The  summation  of  changes  to  the  dependent  variables  within  the 
given  flow  domain  must  then  be  equal  only  to  the  changes  of  the  boundary  nodes 
and  the  overall  contribution  of  any  internal  source  terms. 
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Considering  the  integration  of  the  inviscid  terms  over  the  domain: 


d 

dt 


j  f  UldmdO  +  <f  domain  (Fd9  —  Gdm)  —  f  f  Hidmdd.  (3.57) 

J  J  domain  J  boundary  J  J  domain 


Changes  in  the  dependent  variables  calculated  in  each  cell  are  distributed  equally  to 
each  of  the  corner  nodes  in  a  cell.  Therefore,  the  changes  in  the  nodal  values  over 
the  domain  can  be  expressed  as  a  summation  of  the  changes  within  the  cell  and  the 
boundary  and  source  terms: 


E  W)  =  E 

nodes  cells 


\ 

E  («y 

i  cell  , 

\  nodes  / 


+  boundary  and  source  terms. 


(3.58) 


The  inner  summation  on  the  right  hand  side  is  simply  the  change  in  the  dependent 
variables  within  a  single  cell: 

E  W)  =  A Ui.  (3.59) 

cell 

nodes 

Considering  a  single  cell,  as  in  Fig.  3.10,  a  summation  of  the  boundary  fluxes, 
plus  the  inclusion  of  the  average  value  of  the  source  term  within  the  cell,  comprises 
the  A U{  term: 

At 

A  Ui  =  [fu,  -Fe  +  Fs-  Tn\  +  a  tHi,  (3.60) 

where  T  is  the  flux  through  the  boundary  denoted  by  the  appropriate  subscript. 
For  cells  that  share  a  common  face,  opposite  signs  on  ingoing  and  outgoing  fluxes 
yield  a  cancellation  of  fluxes  along  this  shared  face,  except  for  cells  with  faces  along 
domain  boundaries.  Assuming  that  we  have  a  constant  At /A  ratio  throughout  the 
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Figure  3.10:  Cell  flux  notation. 

domain: 

y:  (SUi)  =  ^  ^  IF  +  XI  HiA  1  +  boundary  terms.  (3.61) 

nodes  domain  nodes  I 

< boundary  / 

Equation  3.61  is  an  equivalent  discrete  representation  of  Eqn.  3.57,  except  for  the 
inclusion  of  the  boundary  and  source  terms.  Therefore,  the  scheme  is  conservative. 
A  similar  proof  may  be  developed  for  the  viscous  terms,  showing  that  they  also 
satisfy  the  conservation  requirement. 


3.7  Solver  Implementation 

Since  unstructured  grids  do  not  use  a  logical  array-like  structure  to  provide  the 
spatial  relationships  between  mesh  points,  a  separate  accounting  method  is  required 
to  define  this  connectivity  pattern.  Unstructured  grid  methods  normally  maintain  a 
list  of  the  nodes  which  comprise  each  individual  cell,  as  well  as  a  connectivity  array 
which  lists  the  cells  surrounding  each  node.  In  the  current  implementation,  this 
usual  construction  is  altered.  Instead  of  applying  the  unstructured  data  model  at 
the  level  of  an  individual  cell,  micro-blocks  of  N  x  N  cells  or  ( N  + 1)  x  (N  + 1)  nodes 
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are  arranged  in  small,  locally  structured  arrays.  These  “patches”  (or  “chunks”  in 
three  dimensions  [22]),  are  arranged  in  a  globally  unstructured  fashion.  By  using 
this  semi-structured  approach,  the  connectivity  array  storage  requirement  is  reduced 
and  a  number  of  additional  advantages  are  found  during  grid  adaptation  (discussed 
in  Chapter  5).  This  method  of  data  decomposition  also  provides  a  natural  format 
for  the  eventual  implementation  of  a  multigrid  convergence  acceleration  technique. 
Successive  levels  of  grid  coarseness  may  be  generated  by  recursively  eliminating 
alternate  grid  lines  in  a  patch.  Extension  of  the  method  to  parallel  processing 
should  be  quite  natural  because  the  solution  domain  is  already  partially  partitioned 
through  the  use  of  the  micro-blocks.  Implementation  of  both  of  these  techniques  is 
recommended  for  future  extensions  of  this  effort. 

The  two  primary  data  structures  used  in  the  code  are  the  CMESH  array,  which 
lists  the  nodes  in  a  particular  patch,  and  the  CNEIB  array,  which  lists  the  neighbor 
patches  surrounding  a  patch.  Globally  maintained  data,  such  as  the  dependent 
variables  and  various  grid  quantities,  are  stored  in  one-dimensional  arrays.  Entries 
in  the  CMESH  array  are  pointers  into  these  variable  arrays;  the  position  of  a  node  in 
the  CMESH  array  provides  the  logical  structure  of  the  patch.  Figure  3.11  illustrates 
the  entries  in  the  CMESH  array.  An  entry  in  a  one-dimensional  variable  array  is 
unique,  implying  that  nodes  on  a  patch  boundary  appear  in  two  CMESH  arrays 
along  edges  and  in  four  CMESH  arrays  in  the  case  of  corner  nodes.  Using  this  data 
structure,  the  calculation  of  fluxes  proceeds  in  a  series  of  three  phases:  data  gather, 
computation,  and  scatter.  Data  used  in  the  flux  calculation  for  a  particular  chunk  are 
first  copied  (gathered)  into  local  patch  arrays,  the  required  computations  performed 
using  these  logically  structured  arrays,  and  the  results  copied  (scattered)  back  into 
the  global  arrays.  This  methodology  insures  that  indirect  array  references,  which  are 
a  serious  performance  penalty,  happen  only  twice,  in  the  gather  and  scatter  steps, 
and  do  not  occur  repeatedly  throughout  the  computation  phase. 
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global  array 


Figure  3.11:  CMESH  array  pointer  entries. 

The  other  data  structure,  the  CNEIB  array,  contains  the  connectivity  pattern 
between  individual  patches;  it  is  here  that  the  “unstructured-ness”  of  the  grid  is 
maintained.  An  entry  in  the  CNEIB  array  for  a  particular  patch  contains  a  list  of 
the  patches  which  share  one  of  the  four  edges  of  the  current  patch.  An  edge  which 
lies  along  a  domain  boundary  receives  an  entry  of  zero  in  the  appropriate  location; 
edges  along  periodic  boundaries  will  contain  the  number  of  the  patch  across  the 
boundary.  The  information  contained  in  this  list  is  only  required  if  adaptation  of 
the  grid  is  to  be  performed.  Knowledge  of  the  neighboring  patches  is  used  to  account 
for  a  difference  in  grid  adaptation  level  along  a  patch  edge  and  to  apply  the  various 
logical  rules  required  during  the  adaptation  process. 

The  solution  domain  boundaries  are  the  only  portions  of  the  grid  where  lists 
of  individual  nodes  are  kept.  Each  boundary  entry  contains  two  elements:  the  first 
is  the  nodal  index  of  the  point  on  the  boundary  itself,  the  second  is  the  index  of  the 
next  point  in  from  the  boundary.  These  inner  points  are  used  for  the  extrapolations 
needed  by  many  of  the  boundary  conditions.  Periodic  boundaries  differ  from  this 
format  in  that  the  two  elements  in  the  boundary  entry  are  the  numbers  of  the  two 
nodes  which  share  the  boundary.  Each  patch  also  includes  information  about  its 
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current  adaptation  level,  the  boundary  condition  type  of  each  of  its  four  edges,  and 
an  indication  of  an  adaptation  level  change  across  the  edge.  Pointers  are  also  kept 
to  indicate  the  parent  and  children  of  a  particular  patch  for  use  in  the  adaptation 
routines. 
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Chapter  4 


TURBULENCE  MODELING 


Two  different  turbulence  models  have  been  implemented  in  the  current  work.  The 
models  are  each  briefly  described  below;  further  details  may  be  found  in  the  appro¬ 
priate  references.  This  description  is  followed  by  a  discussion  of  the  implementation- 
specific  details  of  each  model. 


4.1  Algebraic  Model 

The  algebraic  turbulence  model  employed  is  a  modified  form  of  the  Baldwin- Lomax 
algebraic  eddy  viscosity  model  [7]  (referred  to  hereafter  as  B-L).  Although  this  model 
was  originally  developed  for  steady  flows,  it  has  become  a  de  facto  standard  for  many 
applications  because  of  the  ease  of  its  implementation  and  the  minimal  burden  placed 
on  computational  resources.  The  eddy  viscosity  in  the  boundary  layer  is  modeled 
using  a  two-layer  formulation  which  requires  a  search  normal  to  the  wall  in  order  to 
determine  the  location  of  the  boundary  between  the  layers.  As  is  described  below, 
this  search  causes  a  number  of  complications  in  an  unstructured  grid  implementation. 

In  the  inner  layer  of  the  boundary  layer  model,  the  turbulent  eddy  viscosity  is 
given  by  the  Prandtl-van  Driest  expression: 


VTi  =  p(KiyD)2  |D 


(4.1) 
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where 

i 

-y  { ^  1  —  t)  2  /A+ 

\  Pw  j 

0|  is  the  vorticity  magnitude  which,  for  the  quasi-3D  formulation  is  given  by 

1  ,  0v@  1  0vm  rm 

2  dm  r  06  6  r 

The  variable  y  represents  the  distance  normal  to  the  nearest  blade  surface,  A+ 
is  the  sublayer  thickness,  and  Kq  is  the  von  Karman  constant.  Jobe  and  Hankey 
[33]  have  shown  that  increasing  A+  and  /ci  improves  the  ability  of  the  B-L  model 
to  predict  turbulent  boundary  layers  in  adverse  pressure  gradients.  Optimization 
of  these  parameters  is  beyond  the  scope  of  the  current  effort,  and  the  model  is 
implemented  in  its  standard  form  with  A+  =  26  and  =  0.4. 

In  the  outer  region  of  the  boundary  layer,  the  turbulent  eddy  viscosity  is  given 
by: 

pT0  —  2  Cep  Fmax  Umax  Fkleb-,  (4-4) 

where  Fmax  =  max(y  |  Cl  \  D ),  Fkleb  —  [1  +  5.5(Ca'LEB  y/Vmax)6]  \  Vmax  is  the 
value  of  y  where  Fmax  occurs,  «2  =  0.0168,  Ccp  =  1.6,  and  Ckleb  =  0.3.  The 
turbulence  model  switches  from  the  inner  to  the  outer  model  at  the  first  point  at 
which  A.  P To¬ 
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4.2  One-Equation  Model 

Past  experience  has  shown  that  the  B-L  algebraic  model  overpredicts  the  separa¬ 
tion  region  in  the  high-turning,  high-contraction  blade  passages  of  modern  transonic 
compressor  designs  [60].  An  additional  difficulty  with  the  algebraic  model  is  the 
tracking,  through  any  subsequent  rows,  of  the  turbulent  eddy  viscosity  generated 
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in  upstream  blade  rows.  The  impingement  of  the  upstream  blade  row  wakes  has 
a  considerable  impact  on  the  flow  field  and  performance  of  the  downstream  row, 
particularly  in  the  high-speed  cascade  designs  of  interest.  If  an  algebraic  turbulence 
model  is  used,  as  in  [36,  52,  60],  it  is  unclear  how  the  effects  of  an  upstream  blade 
row’s  eddy  viscosity  would  be  propagated  into  the  downstream  row.  One-  and  two- 
equation  models  are  well  suited  to  this  situation  because  convective  and  dissipative 
terms  are  included  in  the  field  equations,  allowing  the  entire  multi-blade  row  eddy 
viscosity  field  to  evolve  over  time.  Additionally,  these  models  have  a  natural  imple¬ 
mentation  in  an  unstructured  grid  environment;  there  is  no  need  to  perform  a  search 
out  from  the  wall. 

A  number  of  one-  and  two-equation  models  are  enjoying  increased  use  as 
problems  become  more  complex  and  the  availability  of  faster  workstations  helps 
to  ease  computational  limitations.  The  one-equation  models  of  Baldwin  and  Barth 
[8],  Johnson  and  King  [34],  and  Spalart  and  Allmaras  [71],  and  the  two-equation 
models  of  Jones  and  Launder  [35],  Wilcox  [79],  Coakley  [16],  and  Menter  [48]  are  all 
being  used  throughout  the  external  and  internal  aerodynamics  communities.  Each 
of  these  models  have  their  relative  advantages  and  disadvantages.  The  primary 
concerns  in  the  current  effort,  beyond  a  model’s  predictive  abilities,  are  its  com¬ 
putational  efficiency  and  robustness.  During  this  work,  the  models  of  Jones  and 
Launder,  Coakley,  Menter,  and  Spalart  and  Allmaras  have  each  been  implemented. 
However,  the  Spalart- Allmaras  model  has  the  advantage  of  a  reduced  equation  set 
and  has  given  the  least  difficulty  in  its  use. 

The  Spalart-Allmaras  (S-A)  model  has  a  number  of  advantages  that  make  it 
particularly  suitable  for  an  unstructured  grid  implementation.  The  first  is  that  the 
model  is  “local”  in  the  sense  that  the  solution  at  one  point  does  not  depend  on  that 
at  other  points.  This  is  contrasted  with  the  B-L  model  which  requires  a  search  away 
from  the  wall  in  order  to  locate  a  local  maximum.  The  model  does  not  depend  on 
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the  use  of  length  scales,  which  increases  its  generality  and  applicability  to  a  wider 
variety  of  flows  without  case-specific  parameter  tuning.  Finally,  it  has  the  distinct 
advantage  of  not  requiring  grids  which  are  significantly  finer  than  those  required  to 
capture  a  turbulent  velocity  profile.  This  implies  an  ability  to  use  coarser  grids  than 
those  required  by  other  one-  and  two-equation  turbulence  models. 

The  S-A  model  uses  a  single  transport  equation  for  the  turbulent  eddy  viscosity. 
The  turbulent  viscosity  is  given  by: 


vy  =  v 


(4.5) 


where  x  =  v  !ui  and  &  is  the  dependent  variable  of  the  following  transport  equation: 
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In  this  equation, 
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As  before,  |  0  |  is  the  vorticity  magnitude  and  y  is  the  distance  to  the  nearest  wall. 
These  equations  are  those  proposed  in  the  original  S-A  model.  In  their  modeling  of 
dynamic  stall  in  airfoils,  Ko  and  McCroskey  [38]  made  a  number  of  modifications  to 
the  model  in  order  to  improve  convergence  by  preventing  S  from  becoming  locally 
negative.  After  repeated  difficulties  with  instabilities,  the  Ko  and  McCroskey  version 
was  abandoned  in  favor  of  the  original  implementation.  As  suggested  in  [71],  limiters 
have  been  included  on  the  value  of  r  to  insure  that  it  remains  in  the  range  0  —  10. 

Constants  for  the  original  model  given  in  [71]  are  q,i  =  0.1355,  c&2  =  0.622, 
<r  =  2/3,  cw 2  =  0.3,  cw3  =  2,  cvl  -  7.1,  k  =  0.41  and  cvA  -  c&i/k2  +  (1  +  cb2)lcr. 
While  the  original  S-A  model  provides  additional  terms  for  the  prediction  of  the 
transition  location,  these  terms  are  omitted  here;  the  flow  fields  are  considered  to 
be  fully  turbulent. 

The  boundary  conditions  which  were  implemented  in  all  cases  include  a  homo¬ 
geneous  condition  at  the  wall  surface,  the  specification  of  the  freestream  value  at  the 
inlet  plane,  and  an  extrapolation  of  the  viscosity  parameter  at  the  exit.  A  nominal 
freestream  value  of  v  —  ^  is  used,  however,  neither  the  solution  converged  values, 
nor  the  stability  of  the  solution  during  transients,  were  sensitive  to  this  value. 

4.3  Implementation  of  the  Models 

The  following  section  discusses  the  implementation-specific  factors  involved  with 
each  model.  Validation  examples  are  presented  in  Chapter  6. 

4.3.1  Algebraic  Model 

As  previously  described,  the  B-L  algebraic  model  requires  a  search  normal  to  the 
wall  in  order  to  determine  the  boundary  between  the  inner  and  outer  models.  This 
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search  causes  difficulties  in  a  pure  unstructured  implementation  because  of  the  com¬ 
putational  overhead  in  determining  nearest  node  information.  Some  authors  have 
simply  ignored  the  problem  by  using  a  structured  grid  adjacent  to  the  blade  sur¬ 
face  and  including  the  turbulence  model  and  the  turbulent  eddy  viscosity  only  in 
the  inner  grid  [62,  69].  Other  authors  have  overcome  this  difficulty  by  either  over¬ 
laying  a  structured  grid  near  the  blade  surface  [46],  or  by  maintaining  a  search  tree 
data  structure  to  facilitate  the  search  for  the  inner-to-outer  model  transition  loca¬ 
tion  [37].  The  latter  method  is  used  in  the  current  effort.  A  search  line,  called  a 
“thread,”  is  generated  and  stored  for  each  surface  point  after  each  grid  adaptation 
cycle.  Since  the  turbulent  viscosity  is  updated  only  occasionally  (every  50  iterations 
for  steady  flow  cases),  these  threads  are  calculated  and  written  to  a  temporary  disk 
file.  Storing  the  threads  allows  the  search  for  the  inner-to-outer  model  transition  to 
easily  be  performed  along  lines  that  can  be  traced  directly  back  to  the  blade  surface. 
The  turbulent  eddy  viscosity  at  all  those  points  which  are  not  directly  connected  to 
a  blade  surface  (due  to  grid  adaptation),  is  calculated  by  a  recursive  interpolation 
from  points  along  the  threads.  A  drawback  of  the  current  implementation  is  the 
lack  of  a  specific  wake  model,  which  would  track  the  wake  centerline  downstream  of 
the  trailing  edge. 

From  past  experience,  the  B-L  model  tends  to  overpredict  regions  of  separated 
flow.  Experience  has  shown  that  much  of  this  is  due  to  “noise”  in  the  search  away 
from  the  wall  for  the  value  of  Fmax.  A  modification  was  introduced  to  only  include 
values  in  this  search  which  have  a  positive  streamwise  velocity  component.  Ignoring 
regions  of  reverse  flow  causes  the  model  to  yield  larger  values  of  eddy  viscosity  in 
these  areas,  helping  to  suppress  the  formation  of  flow  separation. 


4.3.2  One-Equation  Model 


The  one-equation  model  uses  an  equation  set  which  is  quite  similar  to  that  solved 
for  the  flow  equations,  and  is  therefore  solved  using  an  identical  algorithm.  An  72- 
stage  Runge-Kutta  is  used  to  update  the  turbulent  equations  before  each  iteration 
of  the  flow  equations.  The  code  has  been  implemented  so  that  different  Runge- 
Kutta  schemes  may  be  used  for  the  flow  and  turbulence  equations.  This  general 
method  also  allows  the  use  of  sub-iterations  for  the  turbulence  equations,  which 
may  require  smaller  time  steps  than  the  flow  equations.  In  a  treatment  similar 
to  that  used  for  the  flow  equations,  the  convective  and  source  terms  are  evaluated 
during  each  stage  of  the  Runge-Kutta  step,  while  the  diffusive  terms  are  updated 
during  the  first  stage  and  then  frozen  for  the  remaining  stages.  Second-  and  fourth- 
order  damping  is  included  in  the  Runge-Kutta  scheme,  with  the  coefficients  e 2  and 
e4  of  Eqn.  3.23  taken  directly  from  the  previous  iteration  of  the  flow  equations.  The 
damping  required  is  typically  quite  small.  Extensive  numerical  testing  has  shown 
that  damping  with  coefficients  on  the  order  of  10  —  20%  of  those  used  for  the  flow 
equations  is  adequate  to  suppress  formation  of  any  sawtooth  modes.  For  coarse 
grid  solutions,  a  single  iteration  of  the  same  5-stage  Runge-Kutta  scheme  used  for 
the  flow  equations  is  used.  For  adapted  solutions,  sub-iterations  of  the  turbulence 
model  are  occasionally  required  to  prevent  solution  instabilities.  These  instabilities 
typically  occur  in  highly  loaded  compressor  calculations,  Torming  near  the  trailing 
edge  and  causing  rapid  solution  divergence.  At  most,  two  sub-iterations  of  the  4- 
stage  scheme  are  needed.  Additionally,  a  limiter  is  placed  on  the  viscosity  parameter 
to  insure  that  it  does  not  become  negative;  detected  values  less  than  zero  are  reset 
to  one  half  of  the  freestream  value. 

One  concern  in  the  implementation  of  the  turbulence  model  is  the  calculation  of 
the  wall  distance  function  for  rotor-stator  interaction  cases.  For  an  isolated  cascade 
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Figure  4.1:  Contours  of  wall  distance  function  in  a  closely  coupled  turbine  stage. 

computation,  the  wall  distance  does  not  change  and  can  simply  be  computed  once 
and  saved.  For  a  closely  coupled  stage,  there  may  be  situations  where  the  closest 
surface  point  to  a  particular  node  is  on  a  blade  in  a  different  row  than  that  in  which 
the  node  is  located.  An  example  is  shown  in  Fig.  4.1,  where  point  A  lies  in  the 
downstream  blade  row,  but  it  is  actually  closer  to  the  surface  of  the  upstream  blade 
at  this  particular  moment  in  time.  Although  it  is  possible  to  update  the  wall  distance 
function  as  the  grids  move  past  each  other,  the  turbulence  model  is  sensitive  to  wall 
distance  only  for  those  points  that  fall  within  the  logarithmic  region  of  the  boundary 
layer  [9].  Therefore,  as  in  isolated  airfoil  calculations,  the  wall  distance  is  computed 
at  the  beginning  of  the  computation  and  stored. 
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Chapter  5 


ADAPTATION  TECHNIQUES 


One  of  the  major  benefits  of  unstructured  grid  methods  is  the  ability  to  locally  refine 
the  computational  grid  in  order  to  introduce  additional  grid  points  in  regions  where 
increased  resolution  is  desired.  This  chapter  describes  the  generation  of  the  initial, 
or  base  grid,  the  procedure  for  adapting  the  grid,  and  the  method  in  which  flow 
features  are  detected.  Finally,  the  algorithm  used  to  account  for  the  midface  nodes 
created  at  the  adaptation  interface  is  described. 

5.1  Cascade  Meshes 

All  adaptive  grid  flow  solvers  require  the  generation  of  an  initial  coarse  mesh  from 
which  the  convergence  and  adaptation  process  proceeds.  Unstructured  grid  methods 
that  use  triangular  cells  often  generate  the  initial  grid  using  the  adaptation  algo¬ 
rithm,  with  the  distance  to  solid  surfaces  as  a  weighting  function.  Quadrilateral  cell 
methods  generally  use  structured  grid  generators  to  provide  the  initial  coarse  mesh, 
usually  in  the  form  of  H-grids  for  cascade  calculations. 

A  difficulty  with  using  H-grids,  however,  is  the  poor  resolution  of  leading 
and  trailing  edges  with  significant  curvature.  The  flow  field  through  a  transonic 
compressor  cascade,  and  particularly  the  accurate  determination  of  shock  losses,  is 
strongly  dependent  on  proper  calculation  of  the  strength  of  the  leading  edge  bow 
shock.  In  order  to  provide  better  resolution  of  the  leading  and  trailing  edge  radii  in 
the  transonic  cascades  used  in  the  current  effort,  an  O-grid  is  wrapped  around  the 
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blade  and  blended  into  an  H-grid  in  the  passage.  This  blending  can  be  performed 
in  two  ways.  One  method  is  to  use  a  single  H-grid  in  the  blade  passage,  which 
produces  acceptable  grids,  but  can  cause  excessive  grid  skew  up-  and  downstream  of 
the  blades.  To  improve  the  grid  quality,  additional  H-grid  blocks  can  be  introduced 
up-  and  downstream  of  the  O-grid.  Figure  5.1  shows  a  sample  of  both  of  these  grid 
types  at  the  leading  edge  of  typical  transonic  compressor  cascade.  The  unstructured 
grid  is  generated  using  a  very  flexible,  interactive  graphical  grid  generator  which 
allows  a  great  deal  of  user  control  over  all  grid  parameters.  A  description  of  the  grid 
generator  is  provided  in  Appendix  C.  This  gridding  tool  may  be  used  to  generate 
either  the  simple  or  the  multi-block  form  of  the  blended  0-H  grids. 

Each  of  these  blended  grid  types  create  points  which  are  surrounded  by  either 
five  or  six  cells  rather  than  the  usual  four.  These  special  points  have  a  minor 
effect  on  the  flow  solver  algorithm,  requiring  special  treatment  during  both  the  flux 
integration  and  the  adaptation  process.  This  overhead  is  more  than  compensated 
for  by  the  reduction  in  grid  skewness  near  the  blade  surface,  and  the  ease  in  which 
the  solver  may  be  adapted  to  cascades  with  blunt  leading  edges.  Points  surrounded 
by  more  than  four  cells  require  a  modification  of  the  inviscid  flux  terms  to  account 
for  the  presence  of  the  extra  cells.  During  the  flux  summation,  these  special  nodes 
receive  the  usual  |  distribution  of  fluxes  from  each  of  the  cells  around  the  point.  This 
is  followed  by  a  multiplication  of  the  nodal  flux  value  by  either  |  or  |  depending 
upon  whether  there  are  five  or  six  cells  around  the  node,  respectively. 

5.2  Mesh  Refinement 

The  concept  behind  spatial  grid  adaptation  is  to  adjust  the  local  grid  density  to 
match  the  resolution  required  by  the  flow.  This  involves  increasing  the  grid  density 
in  areas  where  high  resolution  is  required,  such  as  near  shocks  and  in  boundary 
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(a)  Simple  0-H  grid 


(b)  Multi-Block  O-H  grid 


Figure  5.1:  Leading  edge  details  of  blended  0-H  grids. 
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layers,  and  remove  grid  points  from  areas  of  nearly  uniform  flow.  Grid  adaptation 
has  historically  followed  two  paths.  The  first  is  to  retain  the  original  mesh  connec¬ 
tivity  and  only  redistribute  the  existing  grid  points,  "attracting”  them  to  regions  of 
high  flow  gradient  [11,  20].  This  method  has  the  obvious  drawback  that  the  grid 
remains  bounded  by  the  initial  number  of  grid  points.  Additionally,  it  is  often  dif¬ 
ficult  to  sufficiently  manipulate  the  grid  to  track  multiple,  complex  flow  features. 
An  alternate  approach  is  to  locally  embed  additional  cells  in  those  regions  of  the 
flow  requiring  increased  resolution.  This  allows  the  grid  to  completely  match  the 
complexity  of  the  flow  features  with  a  minimum  number  of  points. 

Adaptation  by  refinement  may  be  further  divided  into  two  classes:  grid  regen¬ 
eration  and  cell  division-recombination.  Grid  regeneration,  which  is  the  technique 
often  used  for  triangular  element-based  unstructured  methods,  uses  a  method  such 
as  the  advancing-front  technique  of  Peraire  et  al.  [50],  to  form  a  new  grid  using 
flow  field  features  to  determine  the  local  grid  resolution.  While  this  requires  no 
modification  of  the  solver  algorithm  for  the  adapted  grid,  the  disadvantage  of  this 
method  is  the  time  required  to  regenerate  an  entire  grid.  For  unsteady  solutions,  the 
grid  is  adapted  frequently  and  the  adaptation  method  must  not  produce  significant 
overhead. 

With  the  cell  division  method,  which  is  used  in  the  current  effort,  single  cells 
are  simply  divided  in  one  or  both  directions  depending  on  the  value  of  the  flow 
gradients.  In  triangular  element  meshes,  cell  division  must  often  be  followed  by  a 
local  re- tri angulation  of  surrounding  cells  in  order  to  minimize  the  amount  of  grid 
skew  [9]. 

Cell  division  in  quadrilateral  elements  results  in  the  creation  of  midface  nodes 
at  the  boundary  between  successive  levels  of  grid  refinement.  These  midface,  or 
hanging,  nodes  require  special  treatment  in  the  evaluation  of  the  flux  integrals 
(described  in  Section  5.4  below).  The  adaptation  procedure  also  needs  to  avoid 
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Figure  5.2:  Void  and  island  cells. 


the  generation  of  cell  faces  with  more  than  a  single  midface  node.  This  is  prevented 
by  simply  ordering  the  cells  to  be  divided  from  coarsest  to  finest  adaptation  levels. 

An  advantage  of  the  semi-structured  grid  approach  used  in  the  current  effort, 
over  a  standard  purely  unstructured  grid  method,  is  the  sizable  reduction  in  the  addi¬ 
tional  logic  required  during  adaptation.  In  a  purely  unstructured  implementation, 
such  as  [2,  61],  rules  are  applied  to  eliminate  the  generation  of  any  undivided  cells 
completely  surrounded  by  divided  cells  (“voids”)  or  isolated  divided  cells  (“islands”), 
as  illustrated  in  Fig.  5.2.  While  in  some  cases,  removal  of  these  void  and  island  cells 
is  required  by  the  data  structure  [61],  they  are  generally  undesirable  because  of 
the  overhead  required  by  the  midface  nodes.  In  the  semi-structured  approach  with 
patches  containing  more  than  one  cell,  neither  voids  or  islands  may  exist  in  the  mesh, 
greatly  reducing  this  adaptation  overhead. 

A  standard  method  of  locating  new  grid  points  in  an  adapted  cell  is  to  place 
them  at  the  midpoints  of  the  previous  cell  faces.  Edge  nodes  use  the  data  from  the 
neighboring  corner  nodes,  while  nodes  added  at  a  cell’s  center  use  the  average  of  the 
four  corner  nodes.  Another  advantage  of  the  semi-structured  approach  lies  in  the 
ability  to  maintain  the  grid  stretching  distribution  during  the  division  of  a  patch. 
Within  a  patch,  it  is  possible  to  determine  the  stretching  ratio  of  the  nearby  cells 
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(a)  simple  cell  division  (b)  original  grid  (c)  smooth  cell  division 

Figure  5.3:  Variation  in  stretching  ratio  for  two  forms  of  cell  division. 

and  place  new  node  points  within  the  patch  in  accordance  with  this  pre-existing 
stretching.  This  methodology  allows  the  production  of  refined  cells  with  stretching 
ratios  which  are  half  those  of  the  parent  cells,  thereby  reducing  stretching  error 
as  the  grid  is  refined.  The  improvement  in  stretching  ratios  between  the  current 
method  and  the  use  of  simple  cell  division  is  depicted  in  Fig.  5.3.  The  numbers  in 
each  cell  represent  the  stretching  ratio  between  that  cell  and  the  one  immediately 
below  it.  Grids  used  in  actual  computations  are  often  quite  stretched  around  the 
leading  edge  and  the  use  of  this  alternate  method  of  locating  new  points  produces  a 
rather  dramatic  improvement  in  the  grid  smoothness,  as  seen  in  Fig.  5.4 

Points  on  the  blade  surfaces  are  determined  by  an  interpolating  from  a  geo¬ 
metric  data  set  substantially  finer  than  the  initial  grid.  When  interpolating  surface 
points  in  viscous  grids,  which  have  high  cell  aspect  ratios  near  solid  walls,  regions  of 
high  surface  curvature  are  likely  to  cause  incorrect  placement  of  points  within  the 
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blade  surface,  as  illustrated  in  Fig.  5.5.  For  this  reason,  points  along  added  normal 
grid  lines  are  relocated,  smoothly  moving  them  out  from  the  wall  to  remove  this 
error  condition.  Points  to  be  shifted  are  identified  in  a  manner  similar  to  that  used 
to  create  the  “threads”  used  in  the  algebraic  turbulence  model;  generally  no  more 
than  ten  points  must  be  shifted  to  avoid  the  generation  of  these  incorrect  cells. 

5.3  Feature  Detection 

One  of  the  major  difficulties  to  be  overcome  in  the  use  of  any  adaptive  grid  solver 
is  the  identification  of  the  flow  field  feature  to  which  the  grid  is  to  be  adapted. 
Generally  some  a  priori  knowledge  of  the  flow  field  is  necessary  in  order  to  determine 
which  features  are  of  importance  and  what  function  should  be  used  to  locate  them. 
Static  pressure  gradient  is  typically  useful  for  locating  shocks,  but  will  not  detect 
boundary  or  shear  layers.  Mach  number  and  velocity  gradients  indicate  both  shock 
waves  and  shear  layers,  but  the  large  velocity  gradient  in  shear  layers  often  masks  the 
detection  of  shock  waves,  particularly  in  initial  adaptation  cycles  when  shocks  are 
highly  smeared.  Density  gradient  is  often  a  good  compromise  because  both  shocks 
and  shear  layers  are  detected.  Another  technique,  which  has  been  found  to  be  quite 
useful  in  the  current  effort,  is  flow  feature  detection  based  on  the  gradient  in  the 
computed  turbulent  eddy  viscosity.  This  approach  allows  resolution  of  the  wake 
region  further  downstream  of  the  trailing  edge  than  with  the  use  of  velocity  gradient 
alone.  This  section  discusses  the  detection  of  both  strong  and  weak  features,  and 
considers  the  special  case  of  adaptation  to  unsteady  flows. 

5.3.1  Strong  Features 

One  of  the  primary  reasons  to  adapt  the  computational  grid  is  to  improve  the 
local  grid  resolution  near  primary  flow  features.  Of  the  possible  flow  phenomena 
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in  the  transonic  flows  considered  in  this  thesis,  the  only  true  discontinuities  in  the 
flow  are  shock  waves.  Physical  shocks  have  a  thickness  on  the  order  of  several 
mean  free  paths,  far  beyond  the  resolution  of  any  typical  CFD  code,  and  would 
require  continued  grid  refinement  far  beyond  a  level  which  might  be  considered 
“adequate.”  Additionally,  if  grid  adaptation  were  based  solely  on  the  statistical 
approach  described  below,  shocks  would  tend  to  skew  the  statistics  and  mask  many 
of  the  smoother  features.  For  this  reason,  a  separate  pass  is  made  through  the  solu¬ 
tion  domain  to  detect  those  patches  which  contain  shocks.  Shocks  are  found  by 
marking  those  patches  which  have  a  normalized  second  difference  in  pressure  above 
some  user-specified  threshold.  This  pressure  difference  is  the  same  as  the  pressure 
switch  used  in  the  fourth-order  damping,  give  by  Eqn.  3.24.  Threshold  values  are 
generally  set  at  A P  >  0.02,  although  the  detection  algorithm  is  somewhat  insensi¬ 
tive  to  the  exact  value.  An  example  of  a  transonic  compressor  cascade  pressure  field 
and  the  detected  shock  patches  is  shown  in  Fig.  5.6. 

Another  feature  detection  method  which  has  proven  useful,  detects  those 
patches  along  the  wall  that  have  a  y+  value  at  the  first  point  away  from  the  wall 
below  some  threshold.  For  turbulent  boundary  layers,  which  have  high  surface 
velocity  gradients,  there  is  a  certain  point  beyond  which  continued  resolution  does 
little  to  improve  the  solution.  Adding  more  cells  near  the  wall  serves  only  to  produce 
cells  which  further  limit  the  maximum  time  step  that  may  be  taken  by  the  solution 
algorithm.  Currently,  this  threshold  value  is  set  at  a  y+  =  1.  In  practice,  this  limiter 
does  remove  some  wall  patches,  but  adaptation  often  occurs  in  the  next  layer  of 
patches  away  from  the  wall,  which  usually  contains  the  knee  of  the  velocity  profile. 

5.3.2  Smooth  Features 

Once  the  patches  containing  strong  features  are  removed  from  consideration,  a  sweep 
is  made  through  the  remaining  patches  to  adapt  the  grid  in  regions  of  smooth,  but 
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Figure  5.6:  Transonic  fan  cascade  pressure  field  and  detected  shock  patches. 
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rapid,  flow  variation.  The  method  used  triggers  cell  division  and  recombination  using 
the  undivided  differences  of  tow  user-specified  parameters,  usually  density  and  Mach 
number  or  eddy  viscosity.  Threshold  values  for  cell  division  and  recombination  are 
determined  based  on  the  statistics  of  the  gradient  field,  as  described  by  Aftomis  and 
Kroll  [2], 

In  the  semi-structured  implementation  used  in  the  current  effort,  an  entire 
patch  is  adapted  rather  than  a  single  cell.  Therefore,  the  decision  to  adapt  must  be 
made  by  considering  the  aggregate  flow  gradient  throughout  the  patch.  In  practice, 
once  the  patches  containing  strong  features  are  removed,  a  sweep  is  made  through 
the  remaining  patches  to  determine  the  local  mean  and  standard  deviation  of  the 
trigger  function  in  each  patch.  These  are  combined  to  form  a  global  mean  and 
standard  deviation  and,  from  these,  the  upper  and  lower  threshold  values.  For  a 
patch  to  adapt,  it  is  not  necessary  that  every  cell  within  the  patch  be  tagged  for 
division  or  recombination.  Instead,  some  percentage  of  the  cells  within  the  patch 
must  lie  above  or  below  the  trigger  value.  When  adapting  to  two  quantities,  division 
occurs  when  either  gradient  exceeds  the  divide  threshold;  recombination  only  occurs 
when  the  gradient  for  both  quantities  is  less  than  the  threshold  value.  Local  patch 
threshold  values  are  determined  by  using  the  local  values  of  the  mean  and  standard 
deviation.  For  example,  a  patch  will  divide  if  its  local  threshold  is  above  some 
globally  determined  value: 

T  local  ^  T global  ,  (5.1) 

divide  divide 
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Here,  A /  indicates  the  average  value  of  the  trigger  function  /,  and  o  indicates  its 
standard  deviation.  Ciocai  is  a  factor  determined  from  the  mathematics  of  a  normal 
distribution  curve  and  sets  the  statistical  percentage  of  cells  within  a  patch  that  will 
divide.  Cdivide  is  the  user  specified  threshold  value  based  on  the  global  quantities. 
From  simple  statistics  theory,  the  expressions  for  the  mean  and  standard  deviations 
in  both  the  patch  and  the  global  domain  are  given  by: 
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Generally,  division  is  allowed  when  30%  of  the  cells  within  a  patch  exceed  the 
division  threshold,  yielding  a  Ciocai  ~  0.5.  Similar  expressions  can  be  written  for  cell 
recombination,  requiring  the  specification  of  a  Cfuse  parameter.  Normally,  80%  of 
the  cells  in  a  patch  are  required  to  fall  below  the  fusion  threshold  before  the  patch 
is  allowed  recombine.  Because  of  the  obvious  skew  of  the  distribution  in  a  typical 
adaptation  function  histogram,  as  seen  in  Fig.  5.7,  the  Cdivide  is  usually  set  at  about 
0.8,  while  CfUSe  is  normally  set  at  0.3.  Note  that  the  histogram  in  Fig.  5.7  is  shown 
after  the  smooth  feature  detection,  but  before  the  application  of  the  adaptation 
rules.  Many  of  the  patches  that  satisfy  the  feature  detector  threshold  requirements 
will  not  fuse  or  divide  because  they  violate  these  additional  rules. 
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Figure  5.7:  Typical  adaptation  function  histogram. 

5.3.3  Unsteady  Flow  Detection 

Detection  of  flow  features  in  unsteady  flows  poses  an  additional  challenge  because 
either  the  grid,  the  flow  features,  or  both  are  moving.  The  object  is  to  keep  the 
grid  refined  both  at  the  current  location  of  a  flow  feature,  as  well  as  during  the 
series  of  iterations  before  the  next  adaptation.  Typically,  the  method  used  to  insure 
adequate  resolution  is  to  lower  the  division  threshold  somewhat,  and  to  adapt  the 
grid  frequently.  This  adaptation  must  be  balanced  by  the  increase  in  computational 
time  due  to  the  added  number  of  cells  to  be  processed.  One  method  is  to  adjust  the 
maximum  number  of  adaptation  levels  to  control  generation  of  excessive  cells. 

An  additional  method  is  to  set  a  upper  limit  on  the  number  of  cells  that  can  be 
generated  at  a  given  adaptation  level.  The  cells  allowed  to  adapt  are  ranked  by  their 
adaptation  function  value  and  the  threshold  for  grid  division  adjusted  to  restrain 


84 


Number  of  Patches 


100 

80 

60 

40 

20 

0 

Figure  5.8:  Adaptation  function  histogram  with  cell  division  and  recombination 
thresholds. 

the  continued  growth  of  the  adapted  solution.  Figure  5.8  shows  a  typical  adaptation 
function  histogram  with  the  divide  and  fuse  thresholds,  as  well  as  the  cut-off  value 
used  to  prevent  excessive  cell  production. 

5.4  Numerical  Treatment  of  Midface  Nodes 

The  most  difficult  obstacle  to  the  grid  adaptation  is  the  proper  treatment  of  the  inter¬ 
face  region  between  patches  of  different  adaptation  levels.  The  integration  scheme 
must  maintain  both  conservation  and  accuracy  at  the  interface,  while  not  imposing 
an  undo  computational  burden.  This  section  will  examine  the  problem  in  more  detail 
and  describe  the  algorithm  that  has  been  implemented. 

The  difficulty  at  the  interface  region  is  due  to  two  factors:  a  sudden  change  in 
the  grid  stretching,  or  metric  discontinuity ,  and  the  actual  discontinuity  in  the  grid 
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(a)  Metric  (b)  Grid 

Figure  5.9:  Types  of  grid  interface  discontinuities. 

itself  caused  by  the  midface  nodes  (see  Fig.  5.9).  The  first  of  these  causes  a  significant 
stretching  error,  which  may  harm  the  overall  order  of  accuracy  of  the  method.  On 
the  other  hand,  the  inclusion  of  midface  nodes  requires  special  integration  procedures 
to  insure  that  conservation  is  maintained. 

There  are  typically  two  methods  employed  to  perform  the  flux  integration  at 
the  boundary  between  cells  of  different  adaptation  levels.  The  first  method  includes 
the  hanging  nodes  in  the  integration  of  the  coarse  grid  cells  with  a  custom  integra¬ 
tion  stencil.  Figure  5.10  shows  the  integration  distribution  stencils  for  a  standard 
cell  and  for  cases  of  cells  with  one  or  two  hanging  nodes.  (The  semi-structured 
approach  prevents  any  cells  from  having  more  than  two  hanging  nodes  for  patches 
containing  more  than  one  cell).  This  treatment  has  been  shown  to  be  conserva¬ 
tive,  but  to  suffer  from  an  excessive  stretching  error  [19,  37].  This  interface  method 
also  presents  a  number  of  difficulties  in  its  implementation  in  the  code.  Hanging 
nodes  are  only  contained  within  the  nodal  definition  array  of  the  fine  level  patch 
and  are  not  part  of  the  coarse  patch  data  structure.  Therefore,  inclusion  of  these 
hanging  node  terms  in  the  coarse  cell  integration  would  require  unwanted  commu¬ 
nication  beyond  the  boundary  of  a  single  patch.  A  form  of  this  method  has  been 
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(a)  standard  cell  (b)  one  midface  node  (c)  two  midface  nodes 

Figure  5.10:  Midface  node  cell  integration  stencils  for  conservative  interface 
approach. 

implemented  by  assuming  that  these  hanging  nodes  have  a  value  derived  from  an 
average  of  values  of  the  neighboring  nodes  along  the  edge.  This  requires  a  secondary 
sweep  through  those  coarse  cells  that  border  finer  cells  to  correct  the  fluxes  in  the 
four  corners  of  the  coarse  cell,  without  modifying  the  coarse  cell  flux  contribution 
to  the  hanging  node.  This  method  allows  the  use  of  data  available  within  a  patch, 
without  requiring  communication  to  the  surrounding  patches.  A  final  pass  is  needed 
to  update  all  hanging  nodes  with  the  average  of  their  neighbors.  Use  of  this  con¬ 
servative  interface  approach  is  only  of  value  in  the  evaluation  of  the  inviscid  fluxes, 
which  are  more  susceptible  to  conservation  errors.  Calculation  of  viscous  fluxes  is 
quite  sensitive  to  stretching  error  and  therefore  requires  use  of  an  interface  method 
that  maintains  accuracy  at  the  expense  of  conservation. 

The  operations  in  the  conservative  treatment  may  be  summarized  as  follows: 

•  Complete  a  normal  flux  integration  sweep  through  all  cells  in  the  patch. 

•  Along  each  edge  that  borders  a  finer  cell: 


1.  Calculate  the  midface  node  value  by  averaging  the  neighboring  node 
values. 
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2.  Calculate  the  change  in  the  cell  flux  resulting  from  this  midface  node. 

3.  Using  the  modified  distribution  stencils,  distribute  this  change  to  the 
corners  of  the  coarse  cell. 

®  After  the  global  variable  update,  calculate  the  hanging  node  dependent  vari¬ 
ables  by  averaging  the  values  at  the  adjacent  nodes. 

The  accurate  approach  treats  the  nodes  along  the  interface  as  if  they  belonged 
only  to  the  coarse  cell,  eliminating  the  stretching  error  that  exists  in  going  from  a 
coarse  to  fine  cell.  A  “supercell,”  which  is  a  combination  of  four  fine  grid  cells  as 
showm  by  the  dashed  box  in  Fig.  5.11,  is  used  to  perform  this  updated  integration. 
The  fluxes  calculated  during  this  secondary  pass  replace  the  fine  cell  fluxes  at  the 
coarser  level  cell  nodes  along  the  interface  originally  found  by  considering  the  fine 
cells  only  (i.e.,  only  nodes  A  and  C  in  Fig.  5.11).  This  treatment  ignores  the  presence 
of  the  hanging  nodes  in  the  flux  integration  of  coarse  level  cells  along  the  patch 
boundary.  However,  these  nodes  contribute  to  the  fluxes  in  the  fine  level  patch 
cells.  Once  the  nodal  fluxes  have  been  determined  and  the  dependent  variables  are 
updated,  the  hanging  node  variable  values  are  determined  by  averaging  the  variables 
at  the  adjacent  nodes.  There  are  three  choices  for  the  selection  of  this  averaging 
of  neighboring  nodal  values.  The  first,  which  is  to  perform  a  simple  average  of  the 
dependent  variable  quantities,  must  be  discarded  because  of  the  nonlinearity  of  the 
flow  equations.  The  second  choice,  which  would  allow  the  accurate  treatment  to  also 
be  conservative,  is  to  use  the  dependent  variable  values  resulting  from  an  average  of 
the  flux  vector  values  (F  or  G  in  Eqn.  2.16).  Unfortunately,  this  would  require  the 
solution  of  an  equation  which  is  of  fourth  order  in  the  dependent  variable  values! 
The  final  choice  is  to  average  the  primitive  variables,  (p,vm,v$,p),  which  alleviates 
the  nonlinearity  problem  and  is  computationally  acceptable. 
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Figure  5.11:  Formation  of  “supercell”  for  accurate  interface  treatment. 

The  operations  in  the  accurate  treatment  may  be  summarized  as  follows: 

•  Complete  a  normal  flux  integration  sweep  through  all  cells  in  the  patch. 

•  Along  each  edge  that  borders  a  coarser  cell: 

1.  Re-initialize  the  flux  values  of  the  coarse  nodes. 

2.  Perform  a  flux  integration  in  the  “supercells”  along  this  edge. 

3.  Distribute  the  calculated  fluxes  to  the  coarse  nodes  on  the  edge  only. 

•  After  the  global  variable  update,  calculate  the  dependent  variables  at  hanging 
nodes  by  averaging  the  values  at  their  adjacent  nodes. 

A  similar  operation  is  performed  in  the  calculation  of  the  second-  and  fourth- 
order  smoothing  terms  and  is  also  included  in  the  updates  of  the  one-equation  tur¬ 
bulence  model.  In  addition,  it  was  also  found  necessary  to  perform  the  hanging  node 
averaging  of  the  turbulent  eddy  viscosity. 

A  final  note  in  the  implementation  of  the  accurate  treatment  is  the  impor¬ 
tance  of  detecting  and  treating  inside  corners  between  adaptation  levels,  as  shown 
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Figure  5.12:  Accurate  treatment  of  inside  corner  nodes  taking  into  account  all  sur¬ 
rounding  coarse  level  “supercells.'' 

in  Fig.  5.12.  The  aim  of  the  accurate  treatment  is  to  update  the  coarse  grid  nodes 
along  an  interface  boundary  with  an  integration  stencil  fully  at  the  coarse  grid  level. 
Therefore,  logic  has  been  introduced  into  the  code  to  detect  inside  corners  and  apply 
the  coarse  grid  integration  stencil  to  the  fine  grid  level  patches  that  share  these  corner 
nodes.  When  two  adjacent  edges  of  a  coarse  grid  level  patch  are  found  to  be  divided, 
a  search  is  made  around  the  corner  node  to  tag  all  the  patches  sharing  this  node. 
This  allows  the  algorithm  to  detect  inside  corner  nodes  which  occur  at  the  special 
grid  block  intersection  points  with  five  or  six  surrounding  cells. 
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Chapter  6 


RESULTS 


This  chapter  presents  results  from  the  use  of  the  computational  model  described 
in  the  previous  chapters.  Results  are  first  presented  for  some  of  the  simple  test 
cases  used  to  validate  the  solver.  Computations  of  steady  cascade  flows  in  transonic 
turbomachinerv  are  then  given,  followed  by  an  example  case  using  the  unsteady 
adaptive  solver. 


6.1  Solver  Validation 

A  variety  of  validation  exercises  were  performed  in  order  to  evaluate  various  portions 
of  the  computational  model  separately.  Unless  otherwise  stated,  all  basic  test  cases 
are  performed  with  unit  radius  and  streamtube  thickness.  In  steady  flows,  conver¬ 
gence  is  typically  assumed  once  the  L 2  norm  of  all  dependent  variable  residuals  has 
fallen  four  orders-of-magnitude.  When  using  the  dual  time  stepping  procedure  for 
unsteady  flows,  each  step  in  real  time  is  allowed  to  converge  two  orders-of-magnitude. 
Unsteady  periodic  flows  are  converged  when  an  integrated  quantity,  such  as  lift  or 
drag  coefficient,  has  become  periodic.  Non-periodic  unsteady  flows  are  assumed  to 
be  converged  at  each  real  time  step  as  the  solution  evolves. 
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6.1.1  Viscous  Terms 


The  first  case  is  a  flat  plate  compressible  Blasius  boundary  layer  flow  with  an  inlet 
Mach  number  of  M0 0  =  2.0  and  a  Reynolds  number  based  on  plate  length  of  Rei  = 
104.  This  case  was  computed  with  a  53  x  41  grid  with  an  uniform  axial  spacing,  a 
near-wall  spacing  of  10~3  blade  lengths,  and  a  uniform  geometric  stretching  factor 
of  1.2  in  the  direction  normal  to  the  wall.  Uniform  supersonic  freestream  conditions 
are  imposed  at  the  upstream  boundary,  with  a  no-change  condition  at  the  exit. 
The  upper  surface  uses  a  similiar  radiative  condition,  while  the  wall  surface  uses 
the  standard  adiabatic  no-slip  condition,  with  the  leading  edge  of  the  plate  at  the 
first  grid  point  downstream  of  the  inlet  boundary.  The  number  of  points  in  the 
direction  normal  to  the  wall  was  chosen  to  insure  that  the  Mach  wave  generated  by 
the  wall  leading  edge  singularity  intersects  the  exit  boundary  rather  than  the  upper 
edge  of  the  domain,  as  shown  in  Fig.  6.1.  Figure  6.2  shows  the  velocity  profiles 
and  wall  skin  friction  distribution  along  the  plate.  The  theoretical  curve  is  taken 
from  the  incompressible  Blasius  solution  given  by  White  [78,  pp.  261-267],  with  a 
compressibility  correction  applied  for  the  supersonic  flow  condition  [78,  pp.  586-591]. 
Agreement  wdth  the  classical  Blasius  solution  is  quite  satisfactory. 

6.1.2  Turbulence  Models 

A  flat  plate  in  a  supersonic  flow  with  M*,  =  1.5  and  a  Reynolds  number  based 
on  plate  length  of  Rei,  =  106  was  used  to  validate  the  one-equation  turbulence 
model.  This  case  is  implemented  with  the  supersonic  inlet  and  mixed  exit  conditions 
described  in  Appendix  B.  Cases  were  run  to  show  the  effect  of  spacing  normal  to 
the  wall  on  the  computed  boundary  layer  profile  and  skin  friction.  Three  grids  were 
used,  each  with  53  points  spaced  uniformly  in  the  axial  direction.  The  values  of  the 
first  point  spacing,  normalized  by  the  plate  length,  were  2.5  x  10“%  5.0  x  10~5.  and 
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Figure  6.1:  Blasius  boundary  layer  solution. 


1.0  x  10-4,  which  is  equivalent  to  y+  —  1.77,  2.47,  and  3.54  at  the  exit  of  the  plate.  In 
each  case,  a  uniform  stretching  factor  of  1.2  was  applied  in  the  normal  direction,  and 
the  normal  grid  count  was  adjusted  to  keep  the  upper  boundary  approximately  1.5 
plate  lengths  above  the  wall  surface.  Standard  no-slip  conditions  are  implemented 
along  the  surface  and  a  freestream  condition  along  the  upper  boundary.  As  described 
in  Chapter  4,  a  freestream  value  of  v  —  ^  is  used,  which,  using  Eqn.  4.5,  results  in  an 
imposed  freestream  value  of  v t  =  2.8  x  10~'  v.  Figures  6.3  and  6.4  show  comparisons 
of  the  three  cases  in  terms  of  nondimensional  velocity  profiles  and  skin  friction. 
Both  of  the  theoretical  skin  friction  curves  shown  in  Fig.  6.4  have  been  corrected 
for  the  effects  of  the  =  1.5  flow.  The  laminar  skin  friction  curve  uses  the  same 
correction  as  that  described  in  the  previous  section  for  the  Blasius  boundary  layer. 
The  turbulent  curve  uses  Cj  =  0.455/  In2  0.067^  [78,  pg.  498]  with  the  correction 
of  van  Driest  [73]  found  in  Schlichting  [65,  pg.  717].  These  figures  show  that  the 
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Figure  6.2:  Blasius  boundary  layer  results  for  =  2.0  and  Re^  =  104. 


Skin  Friction  Coefficient 


Figure  6.4:  Skin  friction  distribution  for  three  normal  point  spacings 


S-A  model  gives  results  which  are  quite  acceptable,  even  for  cases  when  there  are 
only  a  few  grid  points  within  the  laminar  sublayer. 

6.1.3  Unsteady  Vortex  Shedding 

The  time  accuracy  of  the  dual  time  stepping  methodology  was  evaluated  by  a  cal¬ 
culation  of  the  well-known  unsteady  vortex  motion  behind  a  circular  cylinder  in  a 
uniform  crossflow.  While  a  substantial  amount  of  effort  has  been  expended  to  com¬ 
pute  this  flow  (see  for  example  Beaudan  [10]),  the  method  described  in  this  thesis 
lacks  the  temporal  or  spatial  accuracy  required  to  capture  the  full  spectrum  of  tem¬ 
poral  and  spatial  scales.  Instead,  the  purpose  of  the  test  was  to  evaluate  the  ability 
of  the  method  to  capture  the  large  scale  motion  through  a  prediction  of  the  frequency 
of  the  vortex  shedding  behind  the  cylinder.  A  row  of  unit  diameter  cylinders  spaced 
twenty  diameters  apart  was  studied  using  an  inlet  Mach  number  of  0.2  and  Reynolds 
number  based  on  diameter  of  1000.  For  these  conditions,  the  Strouhal  number,  as 
given  by  Schlichting  [65,  pp.  31-32],  was  expected  to  be  0.21.  The  multi-block  0- 
H  mesh  used,  shown  in  Fig.  6.5,  contained  a  total  of  4409  nodes  in  265  patches. 
Starting  with  a  uniform  flow  velocity  condition  throughout  the  domain,  the  compu¬ 
tation  tended  to  converge  to  a  steady-state  solution  with  a  stable,  symmetric  vortex 
pattern.  As  the  computation  is  continued,  sufficient  numerical  error  accumulates 
to  make  the  solution  asymmetrical,  causing  alternate  vortex  shedding  to  begin.  In 
order  to  speed  the  formation  of  the  vortex  shedding,  a  surface  blowdng  boundary 
condition  w^as  applied  at  five  points  along  the  upper  surface  of  the  cylinder  for  one 
shedding  period.  This  condition  produces  sufficient  solution  asymmetry  to  trigger 
the  shedding  of  the  first  vortex.  After  the  blowing  is  discontinued,  the  alternate 
vortex  shedding  continues  and  rapidly  becomes  periodic. 

One  of  the  advantages  of  the  dual  time  stepping  method  is  the  ability  to  select 
the  temporal  resolution  of  the  unsteady  solution  through  specification  of  the  real  time 
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Figure  6.5:  Grid  used  for  cylinder  vortex  shedding  test:  detail  of  grid  near  surface 
showing  multi-block  structure. 

step,  At.  This  may  allow  the  prediction  of  large  scale  unsteady  phenomena,  without 
requiring  the  resolution  of  the  smallest  time  scales,  which  are  related  to  the  spatial 
resolution  of  the  grid.  The  drawback,  of  course,  is  that  the  time  step  size  required 
to  yield  an  adequate  level  of  temporal  resolution  is  not  known  a  priori.  To  illustrate 
this  in  the  current  test  case,  the  period  of  the  expected  shedding  was  divided  into  10, 
20,  40,  80,  and  160  equal  real  time  steps.  Because  the  initial  transition  from  uniform 
to  shedding  flow  was  not  of  specific  interest,  the  20,  40,  80,  and  160  time  steps  per 
period  cases  were  started  from  the  partially  developed  10  step  per  period  case  at  a 
nondimensional  time  of  t  =  50,  each  solution  progressing  until  t  =  100,  by  which 
time  all  solutions  had  become  periodic.  At  each  step  in  real  time,  the  solution  was 
converged  in  pseudo-time  until  the  average  of  the  four  dependent  variable  residuals 
had  been  reduced  by  two  orders-of-magnitude. 

The  vortex  shedding  motion  is  quite  evident  in  Fig.  6.6,  which  shows  Mach 
number  contours  for  one  period  of  the  80  iterations  per  period  solution.  The  evo¬ 
lution  of  the  unsteady  lift  and  drag  coefficients  of  the  cylinder  for  the  complete  80 
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iterations  per  period  case  are  shown  in  Figs.  6.7  and  6.8.  Experimental  data  of 
Roshko  [63,  64]  and  Cantwell  and  Coles  [13]  indicate  that  the  drag  coefficient,  for 
Re-D  =  1000  flow,  should  be  Cr>  ~  0.95  —  1.0;  the  computed  results  yield  an  average 
value  of  Cd  ~  1.05,  showing  reasonable  agreement.  A  comparison  of  all  five  iteration 
count  cases  is  presented  in  Figs.  6.9  and  6.10,  which  clearly  showTs  the  considerable 
discrepancy  in  the  computed  shedding  periods  between  the  different  solutions,  an 
effect  of  the  temporal  resolution.  These  figures  also  show  the  asymptotic  behavior  of 
the  solutions;  as  the  iteration  count  is  increased  the  solutions  tend  to  converge  to  a 
final  value.  It  is  also  apparent  from  these  figures  that  there  is  a  minimum  number  of 
iterations  required  to  achieve  even  a  marginally  acceptable  result;  the  10  iterations 
per  period  case  appears  to  be  below  this  allowable  level. 

In  order  to  evaluate  the  effect  of  grid  resolution  on  the  predicted  shedding 
frequency,  solutions  were  also  obtained  with  a  considerably  finer  grid.  This  fine  grid 
was  constructed  by  triggering  adaptation  of  all  patches  in  a  coarse  grid  solution, 
forming  a  computational  domain  with  1060  patches  containing  17297  nodes.  The 
initial  dependent  variable  values  on  the  fine  grid  were  found  by  interpolating  from 
the  same  partially  converged  10  iterations  per  period  case  used  to  start  the  coarse 
grid  cases.  Comparisons  of  the  unsteady  lift  and  drag  coefficients  for  40.  80,  160, 
and  320  iterations  per  period  cases  on  the  fine  grid  are  given  in  Figs.  6.11  and  6.12. 
The  solution  on  the  finer  grid  shows  a  small  increase  in  the  peak  lift  and  a  higher 
overall  drag  coefficient.  A  examination  of  the  velocity  flow  field  downstream  of  the 
cylinder,  as  shown  in  Fig.  6.13,  indicates  the  degree  to  which  smaller  scale  vortices 
are  captured  on  the  finer  grid. 

Summaries  of  the  results  from  all  calculations  on  both  grids  are  provided  in 
Table  6.1,  which  also  includes  the  results  from  the  base  grid  case  using  a  standard 
explicit  time  integration  algorithm,  i.e.  without  the  dual  time  stepping  technique. 
The  explicit  calculation  on  the  refined  grid  was  not  completed  due  to  the  excessive 
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Figure  6.8:  Unsteady  drag  coefficient  for  cylinder  vortex  shedding;  complete  calcu¬ 
lation. 
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Figure  6.9:  Comparison  of  the  unsteady  lift  coefficient  using  different  iteration 
counts;  base  grid. 


Figure  6.10:  Comparison  of  the  unsteady  drag  coefficient  using  different  iteration 
counts;  base  grid. 
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Figure  6.11:  Comparison  of  the  unsteady  lift  coefficient  using  different  iteration 
counts;  refined  grid. 


Figure  6.12:  Comparison  of  the  unsteady  drag  coefficient  using  different  iteration 
counts;  refined  grid. 
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time  requirement:  an  estimated  10  CPU  hours  per  period.  These  results  show  that 
the  computed  Strouhal  number  for  the  coarse  grid  solution  approaches  0.168,  while 
the  fine  grid  values  approach  0.202.  considerably  closer  to  the  expected  experimental 
value  of  0.21.  This  result  seems  quite  reasonable  considering  the  coarseness  of  the 
computational  grid;  additional  levels  of  grid  refinement  would  no  doubt  lead  to 
further  improvement  in  the  computed  Strouhal  number.  It  may  also  be  noted  that 
fewer  pseudo-time  iterations  are  required  to  converge  the  solution  at  each  real  time 
step  as  the  iteration  count  per  period  is  increased.  There  is  less  change  in  the  flow 
field  structure  with  a  smaller  real  time  step,  and  therefore  a  more  rapid  convergence 
to  the  new  state.  Also  included  in  the  tables  is  an  accounting  of  the  effective  CFL 
numbers  used  in  the  implicit  scheme.  These  are  computed  as  the  ratio  of  the  time 
step  in  real  time  to  that  in  pseudo-time,  with  averages  weighted  by  the  cell  areas. 
Finally,  the  CPU  time  requirements  are  provided  for  the  single  processor  Silicon 
Graphics  Power  Indigo2  workstation  used  for  all  the  calculations  in  this  effort.  These 
timing  values  clearly  show  the  considerable  time  savings  possible  with  the  use  of  the 
implicit  method. 
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Base  Grid 


time  steps 
per  period 

iterations 
:  per  period 

computed 
Strouhal  No. 

CFL  Number 
max  avg  min 

CPU  time  per 
period  (min) 

10 

272 

0.144 

7833 

330 

8.0 

2.7 

20 

568 

0.160 

3945 

165 

3.5 

2.6 

40 

800 

0.165 

1961 

82 

1.8 

3.9 

80 

1160 

0.166 

981 

41 

0.9 

5.7 

160 

1646 

0.167 

492 

21 

0.4 

7.3 

expl 

8514 

0.164 

3.5 

203.4 

Refined  Grid 


time  steps 
per  period 

iterations 
per  period 

computed 
Strouhal  No. 

CFL  Number 
max  avg  min 

CPU  time  per 
period  (min) 

40 

2170 

0.205 

8736 

223 

3.5 

47.4 

80 

3402 

0.201 

4390 

111 

1.8 

63.0 

160 

4230 

0.202 

2170 

56 

0.9 

96.6 

320 

5158 

0.202 

1066 

28 

0.4 

122.4 

expl 

25720 

- 

3.5 

614  (est) 

Table  6.1:  Summary  of  circular  cylinder  calculations. 


6.1.4  Intra-Blade  Row  Interface 


The  interface  between  grids  in  relative  motion  was  evaluated  with  a  calculation  of  a 
hot  streak  in  a  uniform  subsonic  flow.  The  base  flow  has  a  uniform  inlet  condition 
with  M<x>  =  0.70  and  Rei  =  103,  where  L  is  the  length  of  each  grid  block  in 
the  streamwise  direction.  The  hot  streak  is  modeled  by  increasing  the  inlet  total 
temperature  by  a  factor  of  1.2  over  a  specified  region.  A  similiar  condition  was  used 
by  Rai  [55]  for  calculation  of  hot  streaks  in  a  turbine  stage.  The  grids  used,  as  shown 
in  Fig.  6.14,  were  both  21  x  21  grids  with  uniform  spacing  in  each  direction.  The 
downstream  grid  was  translated  vertically  upwards  with  a  speed  equivalent  to  half  of 
the  inlet  base  velocity.  Figure  6.15,  showing  temperature  contours  for  a  calculation 
without  grid  motion,  shows  no  evidence  of  the  interface.  This  is  expected  since 
the  interface  points  are  aligned  and  should  not  introduce  any  interpolation  error. 
A  similar  steady  calculation  was  also  performed  with  the  grids  offset  by  half  of 
one  vertical  cell  spacing.  Figure  6.16  shows  the  result  of  calculations  using  two 
and  three  points  for  the  interpolation  at  the  interface.  It  is  evident  that  the  two- 
point  interpolation  introduces  errors  at  the  interface  not  found  when  the  higher- 
order  interpolating  function  is  used.  Similar  calculations  were  performed  with  and 
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Figure  6.15:  Temperature  contours  for  the  coarse  grid  solution  without  grid  motion. 


without  conservation  of  the  interpolated  variables  at  the  interface;  no  detectable 
difference  in  the  solutions  was  found.  A  comment  is  in  order  about  applying  the 
conservative  interface  condition  to  variables  with  an  expected  mean  value  of  zero. 
The  conservative  condition  uses  the  ratio  of  the  integrated  flow  variables,  determined 
before  and  after  the  interpolation,  to  perform  a  correction  of  the  mean  value.  For 
variables  with  a  zero  mean,  this  ratio  should  theoretically  be  a  singularity,  but  in 
practice  will  be  a  ratio  of  two  very  small  quantities.  These  numbers  are  both  on  the 
order  of  the  roundoff  error,  and  their  ratio  can  therefore  vary  wildly.  Logic  has  been 
introduced  to  detect  this  condition  and  add  an  arbitrary  constant  to  both  integrated 
values  in  order  to  avoid  this  singularity. 

Unsteady  calculations  were  performed  on  the  same  case  in  order  to  validate  the 
rotational  terms  in  the  flow  equations,  as  well  as  the  performance  of  the  interface 
in  unsteady  flows.  As  described  in  the  preceding  section,  the  required  temporal 
resolution  was  determined  by  performing  the  calculation  with  a  varying  number  of 
real  time  steps  per  grid  passing  period.  The  temperatures  at  points  in  the  center  of 
the  downstream  grid  (denoted  by  A  and  B  in  Fig.  6.14)  was  monitored  to  determine 
periodicity.  Using  a  time  step  equivalent  to  20  iterations  per  period,  the  calculation 
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(a)  2  point  interpolation 


(b)  3  point  interpolation 


Figure  6.16:  Temperature  contours  for  solutions  with  different  interface  interpolation 
functions. 
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Figure  6.17:  Temperature  history  for  solutions  with  different  iteration  counts  per 
grid  passing  period. 

was  first  converged  to  a  steady  state  solution  without  grid  motion.  Calculations 
were  then  run  an  additional  six  grid  passings  with  20,  40,  and  80  iterations  per 
passing  period.  The  20  iteration  case  shows  a  discrepancy  between  the  steady  state 
temperatures  and  the  peak  temperatures  in  the  unsteady  case;  this  behavior  is  absent 
from  the  40  and  80  iteration  cases.  Figure  6.18  shows  the  temperature  contours  at 
the  end  of  the  six  grid  cycles.  The  20  iteration  case  shows  a  distortion  of  the  hot 
streak  in  the  moving  downstream  grid,  which  is  improved  as  the  iteration  count  per 
period  is  increased.  The  steadiness  of  the  profile  is  particularly  good  considering  the 
coarseness  of  the  grids. 

6.1.5  Adaptation  to  Unsteady  Flows 

A  test  of  the  unsteady  adaptation  was  performed  for  the  hot  streak  case  described 
above.  At  each  step  in  real  time,  the  grid  was  adapted  to  undivided  differences 
in  temperature  after  the  average  L2  norm  of  the  solution  residuals  had  fallen  a. 
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(c)  80  iterations  /  period 


Figure  6.18:  Comparison  of  temperature  contours  for  coarse  grid  solutions  with 
different  iteration  counts  per  grid  passing  period. 
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single  order  of  magnitude.  The  solution  was  then  allowed  to  converge  another  three 
orders-of-magnitude,  based  on  the  solution  residual  immediately  following  the  adap¬ 
tation.  In  general,  four  iterations  in  pseudo-time  were  performed  before  the  grid  was 
adapted,  followed  by  another  seven  to  complete  the  real  time  iteration  step.  Control 
over  the  total  number  of  grid  points  was  maintained  by  allowing  only  two  levels  of 
cell  refinement  beyond  the  base  grid.  An  additional  solution  using  a  globally  fine 
grid  was  also  obtained  in  order  to  compare  to  the  adapted  solution.  Although  the 
number  of  points  changed  at  each  iteration,  the  adapted  grids  typically  contained 
approximately  460  patches  and  7800  nodes  as  opposed  to  the  fine  grid  of  800  patches 
and  over  13000  nodes. 

There  are  a  number  of  differences  to  note  between  the  original  coarse  grid 
solution  and  that  obtained  using  the  medium,  fine,  and  adapted  grids.  The  first  is 
that  the  enhanced  grid  resolution  lessens  the  numerical  dissipation  and  consequently 
reduces  the  spreading  of  the  hot  streak.  This  is  illustrated  in  Fig.  6.19  which  shows 
the  comparison  of  the  coarse,  medium,  and  fine  grid  solutions.  This  improved  res¬ 
olution  is  also  reflected  in  the  change  in  the  peak  temperature  shown  in  Fig.  6.20, 
which  includes  an  overlay  of  the  fine  grid  steady  solution  values.  Figure  6.20  also 
shows  that  the  adapted  solution  includes  a  small  dip  or  “trough”  in  the  tempera¬ 
ture  distribution  on  one  side  of  the  streak.  An  attempt  was  made  to  remove  this 
anomaly  by  increasing  the  number  of  iterations  per  grid  passing.  Figure  6.21  shows 
a  close-up  of  the  temperature  distribution  for  40,  80,  and  160  iterations  per  grid 
passing  overlaid  with  the  fine  grid  values.  Increasing  the  number  of  time  steps  does 
reduce  the  size  of  the  dip,  but  it  is  equally  possible  that  this  is  an  artifact  of  the 
adaptation  itself;  this  phenomenon  currently  remains  unexplained. 

Figure  6.23  shows  a  selected  number  of  different  images  of  the  grid  at  various 
times  in  a  grid  passing  period.  Improvement  in  the  resolution  of  the  streak  is  evident 
over  that  shown  in  the  previous  images  without  the  adaptation.  Slight  distortion  of 
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the  contours  is  seen  in  the  moving  grid  region,  and  is  likely  due  to  “noise”  in  the 
location  of  the  edge  of  the  region  of  finest  cells.  Some  improvement  in  the  solution 
can  be  made  by  reducing  the  cell  division  threshold,  but  this  results  in  the  finest 
level  of  adaptation  in  all  patches  in  the  downstream  grid,  defeating  the  purpose 
of  the  adaptive  algorithm.  There  is  also  some  evidence  of  the  error  introduced  by 
the  hanging  node  interfaces,  particularly  where  an  internal  or  external  corner  exists 
between  adaptation  levels.  These  corners  cause  an  accumulation  of  stretching  error 
in  both  directions,  resulting  in  a  “worst  case”  test  for  the  interface  algorithm. 
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;ure  6.20:  Temperature  history  for  fine,  medium,  and  adapted  grid  solutions 


Nondimensional  Time 


Figure  6.21:  Comparison  of  unsteady  adaptation  solutions. 
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Figure  6.22:  Comparison  of  adapted  grids  and  temperature  contours  for  three  con¬ 
secutive  periods. 
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6.1.6  Isolated  Cascade  Cases 


Two  single  blade  row  cases  are  presented  in  order  to  validate  the  computational 
model  against  available  experimental  data.  Each  of  these  cases  are  representative  of 
typical  designs  found  in  modern  transonic  turbomachinerv. 

Transonic  Fan  Cascade 

The  compressor  cascade  used  is  a  near-tip  section  from  a  transonic  fan  design  with 
experimental  data  for  the  surface  pressure  distribution  available  in  Fleeter  et  ah  [25]. 
The  case  chosen  has  an  inlet  relative  Mach  number  of  1.535  and  a  Reynolds  number 
based  on  blade  chord  and  inlet  relative  conditions  of  1.16  x  106:  experimental  data 
are  available  for  a  static  pressure  ratio  range  of  1.22  -  2.30. 

Calculations  were  performed  for  a  back  pressure  ratio  of  1.505  with  the  experi¬ 
mentally  determined  streamtube  contraction  (defined  as  the  ratio  of  the  exit-to-inlet 
streamtube  thickness)  of  approximately  96%.  and  with  values  4%  above  and  below 
this  value;  the  streamtube  radius  was  held  const  ant.  In  each  case  the  computations 
proceed  on  an  initial  grid  until  the  positions  of  the  bow  and  passage  shocks  have 
stabilized.  With  the  low  resolution  of  the  typical  initial  grid,  the  passage  shock  can 
be  smeared  over  approximately  20%  of  the  blade  surface.  The  initial  adaptation 
cycles  are  performed  with  fairly  low  gradient  thresholds,  serving  to  stabilize  the 
shock  position  and  reduce  the  smearing  of  the  normal  shock.  The  change  in  shock 
position  and  definition  after  two  adaptation  cycles  may  be  seen  by  comparing  the 
original  and  the  adapted  grids  and  flow  fields  given  in  Fig.  6.24. 

A  comparison  of  the  three  computed  blade  surface  pressure  distributions  to 
that  measured  by  Fleeter,  et  al.  [25]  is  seen  in  Fig  6.25.  Good  agreement  with 
experiment  is  obtained  at  the  design  contraction  ratio,  and  a  reduction  in  the  pas¬ 
sage  shock  strength  with  decreased  contraction  ratio  is  found  as  expected.  It  is 
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Contraction  Ratio 


Suction  Side  X/C  Pressure  Side 


Fi  gure  6.25:  Comparison  of  surface  pressure  distribution  for  three  streamtube  con¬ 
traction  ratios. 

interesting  to  note  that  since  the  overall  cascade  pressure  ratio  is  fixed,  the  pressure 
rise  is  redistributed  within  the  cascade  depending  upon  the  specified  contraction. 
Therefore  a  cascade  with  a  lower  contraction  ratio  will  have  a  stronger  passage 
shock  and  a  weaker  leading  edge  bow  shock.  Comparison  of  the  total  pressure  ratio 
distribution  in  the  exit  wake  at  a  traverse  location  25%  chord  downstream  of  the 
trailing  edge  is  given  in  Fig.  6.26.  Best  agreement  with  experiment  is  found  for  the 
design  contraction  ratio  both  for  wake  shape  and  for  the  mass  averaged  losses.  The 
influence  of  contraction  ratio  on  the  passage  shock  strength  is  also  seen  here  in  the 
shift  in  position  of  the  wake  centerline  indicating  the  variation  in  turning  between 
shocks  of  different  strengths. 
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Figure  6.26:  Experimental  and  computed  total  pressure  ratio  25%  chord  downstream 
of  trailing  edge. 

Transonic  Turbine  Rotor 

Computations  are  also  compared  with  experimental  results  for  the  transonic  turbine 
rotor  cascade  of  Ashworth  et  ah  [6].  The  comparison  case  is  the  nominal  design 
point,  with  an  inlet  Mach  number  of  0.59,  an  incidence  angle  of  58.06  degrees,  a 
Reynolds  number  of  9.19  x  10°  based  on  the  blade  chord,  and  an  isentropic  exit  Mach 
number  of  1.18.  The  calculations  include  the  test  cascade  expansion  of  approximately 
12%  and  use  an  isothermal  wall  specification  with  To/Twau  =  1.5.  The  adaptation 
cycle  used  is  similiar  to  that  described  for  the  compressor  flows  above.  The  initial 
and  adapted  grids  and  flow  fields  are  presented  in  Fig.  6.27,  clearly  showing  the 
improvement  in  the  prediction  of  the  blade  boundary  layer  and  wake  over  that  in 
the  initial  grid  solution.  The  comparison  with  the  experimentally  measured  blade 
loading  in  Fig.  6.28  shows  a  reasonable  agreement  with  the  available  data.  The 
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experimental  and  computed  heat  transfer  distributions  are  given  in  Fig.  6.29,  with 
the  Nusselt  number  defined  as  [6]: 


(JwCq 

(  T0  -  Tw  )k 


(6.1) 


where  qw  is  the  heat  transfer  at  the  wall  and  eg  is  the  tangential  chord.  Figure  6.29 
highlights  the  improvement  in  the  predicted  values  as  the  resolution  of  the  blade 
boundary  layer  is  improved. 
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Figure  6.28:  Turbine  rotor  blade  loading  diagram  for  the  adapted  solution. 


TE  LE  TE 

Figure  6.29:  Turbine  rotor  heat  transfer  distribution  for  the  adapted  solution. 
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6.2  Turbine  Stage  Application 


In  order  to  demonstrate  the  complete  computational  model,  calculations  of  the 
unsteady  interaction  in  a  transonic  turbine  stage  are  presented.  This  turbine  is  a 
single  stage  design  with  36  guide  vanes  and  61  rotors  blades:  the  rotor  is  the  same  as 
that  presented  in  the  previous  section.  Design  parameters  for  the  full  scale  machine 
include  a  pressure  ratio  of  4.3  across  the  stage,  a  vane  isentropic  exit  Mach  number 
of  1.1S,  and  a  wall/gas  temperature  ratio  of  0.63.  The  Reynolds  number  based  on 
the  guide  vane  true  chord  and  isentropic  vane  exit  conditions  is  Re  —  2.47  x  106. 
The  actual  machine  uses  film  cooling  in  both  the  vane  and  rotor  rows  which  is  not 
included  here.  The  section  modeled  is  at  a  near  mid-span  location  and  includes 
an  assumed  radius  and  streamtube  thickness  change  taken  from  the  stage  endwall 
geometry  and  previous  calculations  on  the  stage  [1].  Although  the  computer  code 
described  here  has  been  written  in  a  general  form,  allowing  any  number  of  blade 
rows  to  be  calculated,  each  with  an  integer  number  of  passages,  a  1:1  blade  ratio 
grid  has  been  used  to  reduce  the  overall  computational  requirements.  To  maintain 
the  original  blade  solidity  (ratio  of  blade  chord  to  cirumferential  spacing),  the  down¬ 
stream  blade  row  has  been  uniformly  stretched  by  the  ratio  of  the  row  blade  counts 
in  the  actual  machine  (61/36).  All  computations  use  the  axial  spacing  between  blade 
rows  from  the  original  design.  The  1:1  blade  ratio  grid  used  in  the  following  series  of 
calculations  is  shown  in  Fig.  6.30,  which  includes  a  detail  of  the  gap  region  between 
the  blades. 

The  unsteady  lift  and  drag  coefficients  are  again  used  to  determine  when  con¬ 
vergence  to  a  periodic  solution  has  been  attained.  The  lift  and  drag  forces  are 
resolved  in  the  circumferential  and  axial  directions  respectively,  not  in  a  coordinate 
system  aligned  with  blade  stagger  angle.  The  evolution  of  the  unsteady  vane  and 
rotor  lift  and  drag  coefficients  during  a  complete  calculation  is  shown  in  Figs.  6.31 


for  both  lamina r  and  turbulent  cases,  the  turbulent  solution  evolving  from  a  par¬ 
tially  converged  laminar  solution.  Both  coefficients  are  defined  in  the  standard  form, 
nondimensionalizing  the  integrated  forces  by  the  inlet  vane  upstream  density  and 
velocity  and  the  respective  blade  chords.  The  large  coefficient  magnitudes  are  a 
result  of  the  low  inlet  velocity  (M co  =  0.14)  and  the  considerable  blade  forces  on 
a  machine  with  a  pressure  ratio  of  4.3  and  more  than  70°  of  turning  in  each  row. 
This  choice  of  nondimensionalization  also  leads  to  the  negative  lift  coefficient  on 
the  rotor,  which  is  moving  downwards.  Figure  6.32  provides  a  method  to  reference 
variations  in  the  lift  coefficient  to  the  blade  row  motion:  Mach  number  contours  for 
the  same  four  blade  positions  are  presented  in  Fig.  6.33.  From  these  figures  the 
primary  effects  of  the  blade  motion  on  the  blade  forces  may  be  described.  At  the 
beginning  and  end  of  the  blade  passing  cycle,  (points  A  and  E  in  Figures  6.32  and 
6.33),  the  rotor  passage  is  aligned  with  the  vane  passage,  causing  the  least  blockage 
to  the  flow.  Because  the  stage  mass  flow  is  essentially  constant  throughout  the  cycle, 
the  low  blockage  causes  the  flow  velocity  along  the  blade  suction  surfaces  to  drop, 
lowering  the  lift  forces.  Similarly,  at  position  C,  both  the  rotor  alignment  with  the 
center  of  the  vane  passage  and  the  vane  wake  flow  combine  to  maximize  the  blockage. 
This  increases  the  flow  velocity  across  the  suction  sides  of  both  blades,  raising  the 
lift  forces.  While  the  vane  lift  force  is  primarily  affected  by  the  passing  potential 
field  of  the  rotor,  the  vane  wake  produces  an  additional  effect  on  the  rotor  lift.  This 
is  evident  at  position  D,  where  the  vane  wake  fluid  forms  a  large  vortex  along  the 
rotor  pressure  surface,  further  increasing  the  passage  blockage. 

A  number  of  issues  must  be  resolved  when  calculating  a  realistic  geometry, 
particularly  considering  the  effect  on  the  required  computational  time.  The  object 
of  any  calculation  is  to  provide  the  most  accurate  representation  of  the  physical 
system  in  the  least  computational  time.  Therefore,  it  is  important  to  understand 
the  sensitivity  of  the  computed  solution  to  the  choice  of  the  various  computational 


(a)  lift  coefficient 


(b)  drag  coefficient 


Figure  6.31:  Evolution  of  lift  and  drag  coefficients  for  fully  turbulent  solution 
(50  iterations  per  passing,  2  orders-of-magnitude  residual  drop). 
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parameter.  As  will  be  shown  below,  each  of  these  choices  has  an  effect  on  both  the 
solution  and  the  time  required  for  the  computation.  These  issues  will  be  discussed 
using  calculations  with  the  coarse  grid  shown  in  Fig.  6.30,  before  adapted  grid  results 
are  presented.  The  questions  to  be  considered  include  the  following: 

@  How  many  real  time  steps  are  required  per  blade  passing  period? 

®  How  well  do  the  results  from  the  implicit  solution  compare  to  a  fully  explicit 
solution? 

@  How  many  orders-of-magnitude  drop  in  the  residual  are  required  at  each  real 
time  step? 

•  How  does  the  implementation  of  the  turbulence  model  affect  the  solution? 

The  first  set  of  calculations  were  performed  to  examine  the  effect  of  the  tem¬ 
poral  resolution  on  the  computed  solution.  Coarse  grid  solutions  were  obtained 
for  25,  50.  100,  200,  and  400  real  time  iterations  per  rotor  passing,  each  solution 
obtained  with  two  orders-of-magnitude  residual  drop  per  iteration.  Both  laminar 
and  turbulent  calculations  were  performed,  with  all  solutions  starting  from  a  nearly 
converged  50  iteration  per  passing  solution  and  allowed  to  run  until  periodicity  was 
attained.  A  comparison  of  the  unsteady  lift  and  drag  coefficients  on  the  vane  and 
rotor  is  shown  in  Figs.  6.34  and  6.35  for  the  laminar  calculation,  and  in  Figs.  6.36 
and  6.37  for  turbulent  flow.  As  noted  in  the  unsteady  vortex  shedding  and  hot 
streak  calculations  presented  in  the  previous  sections,  further  reduction  in  the  real 
time  step  size  beyond  200  iterations  per  passing  has  little  effect  on  the  solution. 

In  order  to  evaluate  the  accuracy  of  the  implicit  solver,  fully  explicit  calcula¬ 
tions  were  also  performed.  With  the  density  of  the  grid  used  in  these  calculations, 
the  global  minimum  time  step  was  approximately  850  iterations  per  blade  passing. 
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Figure  6.34:  Variation  of  vane  and  rotor  lift  coefficient  for  different  iteration  counts 
laminar  solution. 
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Figure  6.35:  Variation  of  vane  and  rotor  drag  coefficient  for  different  iteration  count 
laminar  solution. 
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Figure  6.36:  Variation  of  vane  and  rotor  lift  coefficient  for  different  iteration  counts; 
turbulent  solution. 
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Figure  6.37:  Variation  of  vane  and  rotor  drag  coefficient  for  different  iteration  counts 
turbulent  solution. 


An  additional  explicit  calculation,  with  2500  iterations  per  passing  was  also  per¬ 
formed  to  investigate  the  level  of  convergence  of  the  explicit  solution.  A  comparison 
of  the  explicit  solutions  with  the  results  of  the  50  and  200  iteration  per  passing  cases, 
for  both  lift  and  drag  coefficients,  is  presented  in  Figs.  6.38  and  6.39.  These  results 
indicate  a  quite  acceptable  agreement  between  the  implicit  and  explicit  formulations, 
with  little  change  in  the  explicit  solution  as  the  time  step  count  is  increased.  A  com¬ 
ment  about  the  temporal  accuracy  of  the  solution  is  appropriate  here.  The  modified 
multi-stage  Runge-Kutta  algorithm  used  in  the  explicit  algorithm  is  of  first,  order 
temporal  accuracy,  while  the  implicit  formulation  uses  a  second-order  backwards  dif¬ 
ference  temporal  operator.  Therefore,  an  explicit  solution  with  2500  iterations  per 
passing  has  the  same  temporal  accuracy  as  the  implicit  solution  with  50  iterations 
per  passing.  Comparing  solution  CPU  times  for  computations  with  the  same  level  of 
temporal  accuracy ,  there  is  roughly  a  factor  of  8  savings  in  computer  time  between 
the  implicit  and  explicit  solutions. 

The  statistics  of  all  of  these  computations,  shown  in  Table  6.2,  reflect  similar 
trends  as  those  found  during  the  computcition  of  the  cylinder  vortex  shedding  cases. 
The  number  of  pseudo-time  iterations  per  real  time  step  falls  as  the  real  time  incre¬ 
ment  is  reduced.  It  is  also  interesting  to  note  that  the  turbulence  model  has  a  small 
effect  on  the  number  of  pseudo-time  iterations  required  to  converge  the  solution 
at  each  real  time  step,  the  increased  damping  introduced  by  the  turbulence  model 
tending  to  hasten  convergence. 

The  next  issue  relates  directly  to  the  computer  time  required  to  obtain  a.  useful 
solution.  If  the  residual  drop  used  for  each  real  time  step  does  not  greatly  affect  the 
solution,  then  considerably  fewer  pseudo-time  iterations,  and  therefore  less  compu¬ 
tational  time,  will  be  required  for  each  blade  passing  period.  In  order  to  evaluate 
this  choice,  four  solutions  were  obtained,  each  with  200  real  time  iterations  per  blade 
passing,  at  0.5,  1,  1.5,  2,  3,  and  4  orders-of-magnitude  residual  drop.  A  comparison 
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Figure  6.38:  Comparison  of  lift  coefficient  for  implicit  and  explicit  solutions:  laminar 
solution. 
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Figure  6.39:  Comparison  of  drag  coefficient  for  implicit  and  explicit  solutions;  lam¬ 
inar  solution. 


Laminar  Solutions 

real  time 

pseudo-time 

average 

CFL 

CPU  time 

iterations 

iterations 

iterations 

Number 

per  period 

per  period 

per  period 

per  time  step 

max 

avg 

min 

(min) 

25 

217 

8.68 

232 

27.7 

1.86 

1.5 

50 

316 

6.33 

116 

13.8 

0.93 

2.4 

100 

415 

4.15 

58 

6.9 

0.46 

3.2 

200 

547 

2.73 

29 

3.5 

0.23 

4.5 

400 

821 

2.05 

14 

1.7 

0.12 

6.2 

expl 

850 

1.00 

3.5 

6.2 

expl 

2500 

1.00 

3.5 

18.2 

Turbulent  Solutions 

real  time 

pseudo-time 

average 

CFL 

CP  LI  time 

iterations 

iterations 

iterations 

Number 

per  period 

per  period 

per  period 

per  time  step 

max 

a.vg 

min 

(min) 

25 

240 

9.59 

259 

30.4 

1.86 

2.3 

50 

328 

6.56 

130 

15.2 

0.92 

3.0 

100 

405 

4.05 

65 

7.6 

0.46 

3.9 

200 

513 

2.57 

32 

3.8 

0.23 

5.6 

400 

803 

2.01 

16 

1.9 

0.12 

8.8 

expl 

850 

1.00 

3.5 

9.3 

expl 

2500 

1.00 

3.5 

27.4 

Table  6.2:  Results  of  iteration  count  per  period  tests:  2  orders-of-magnitude  residual 
drop. 


of  the  lift  and  drag  coefficients  for  the  0.5.  2,  and  4  orders-of-magnitude  cases  are 
shown  in  Fig.  6.40  and  6.41.  These  results  clearly  show  that  the  choice  of  residual 
does  not  have  a  great  influence  on  the  eventual  solution.  Additionally,  no  influ¬ 
ence  of  the  residual  drop  level  was  observed  on  the  number  of  blade  passing  cycles 
required  to  achieve  a  periodic  solution;  all  solutions  were  essentially  periodic  after 
approximately  20  periods.  The  increased  effort  required  to  attain  the  additional  con¬ 
vergence  is  reflected  in  the  difference  in  the  iteration  count  per  period  for  the  various 
cases  shown  in  Table  6.3.  This  table  also  includes  an  accounting  of  the  relative  error 
magnitude  averaged  over  a  complete  period.  This  error  is  defined  as  the  difference 
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in  force  coefficient  values  between  a  solution  of  some  level  of  convergence  and  the 
solution  with  four  orders-of-magnitude  residual.  These  results  led  to  the  selection 
of  two  orders-of-magnitude  residual  drop  for  many  of  the  other  calculations  in  this 
section. 


orders  of 
residual 
drop 

pseudo-time 
iterations 
per  period 

average 

iterations 
per  time  step 

average  relative  error  (%) 

CPU  time 
per  period 
(min) 

vane 

rotor 

Cl 

Cd 

CL 

CD 

0.5 

400.3 

2.00 

0.029 

0.041 

0.120 

0.132 

4.79 

1.0 

400.9 

2.01 

0.029 

0.041 

0.120 

0.131 

4.80 

1.5 

406.8 

2.03 

0.028 

0.041 

0.119 

0.130 

4.87 

2.0 

513.0 

2.57 

0.017 

0.024 

0.082 

0.126 

5.64 

3.0 

765.7 

3.83 

0.005 

0.007 

0.024 

0.023 

7.40 

4.0 

1249.5 

6.25 

0.0 

0.0 

0.0 

0.0 

11.10 

Table  6.3:  Results  of  magnitude  of  residual  drop  tests;  turbulent  solutions. 


The  final  issue  to  be  examined,  is  the  effect  of  the  turbulence  model  imple¬ 
mentation  on  the  solution.  For  this  example,  laminar  and  turbulent  solutions  are 
compared,  with  an  additional  solution  obtained  using  the  one-equation  turbulence 
model,  but  without  propagation  of  the  turbulent  eddy  viscosity  across  the  interface 
between  blade  rows.  The  no-propagation  case  uses  a  standard  no-change  boundary 
condition  on  the  eddy  viscosity  at  the  exit  of  the  upstream  row,  with  the  eddy 
viscosity  reinitialized  at  the  inlet  of  the  downstream  row.  All  calculations  are  per¬ 
formed  with  200  iterations  per  blade  passing  and  two  orders-of-magnitude  residual 
drop.  The  results  of  this  comparison,  given  in  Figs.  6.42  and  6.43.  show  that,  as 
expected,  the  vast  majority  of  the  difference  between  the  laminar  and  turbulent 
solutions  may  be  attributed  to  the  modeling  of  the  turbulent  boundary  layers  along 
the  individual  blades.  The  propagation  of  the  turbulent  eddy  viscosity  from  the 
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Figure  6.40:  Effect  of  residual  drop  on  turbine  stage  lift  coefficient  (200  iterations 
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6.41:  Effect  of  residual  drop  on  turbine  stage  drag  coefficient  (200 
ssing). 


upstream  to  the  downstream  row  has  only  a  minor  effect  on  the  vane,  with  a  more 
significant  effect  on  the  rotor. 

The  final  figures  in  this  section  are  provided  to  show  the  relative  level  of 
unsteadiness  in  the  solution  and  its  impact  on  blade  surface  parameters.  Figures  6.44 
and  6.45  show  the  blade  loading  and  Nusselt  number  distributions  on  both  blades 
at  ten  equally  spaced  times  during  a  single  blade  passing  period.  As  expected,  the 
unsteadiness  on  the  vane  is  confined  to  the  uncovered  portion  of  the  blade  along 
the  suction  surface,  upstream  of  the  trailing  edge.  Little  effect  is  seen  on  the  pres¬ 
sure  surface,  which  is  not  exposed  to  the  potential  field  of  the  passing  rotor.  The 
distributions  on  the  rotor  clearly  show  the  progression  of  the  vane  wake  as  a  wave 
which  progresses  along  the  blade  surface.  The  wake  motion  is  also  apparent  in  the 
uncovered  portion  of  the  rotor  suction  surface  just  downstream  of  the  leading  edge. 
These  figures  highlight  the  importance  of  considering  not  only  the  time-averaged 
values,  but  also  the  entire  envelope  of  the  unsteadiness,  particularly  for  the  heat 
transfer  prediction. 
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Figure  6.43:  Comparison  of  drag  coefficient  for  laminar,  turbulent,  and  turbulent 
without  eddy  viscosity  propagation  solutions:  200  iterations  per  passing. 


Figure  6.44:  Unsteady  blade  loading  diagrams  throughout  a  rotor  passing  period. 
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Figure  6.45:  Unsteady  blade  heat  transfer  throughout  a  rotor  passing  period. 


6.2.1  Adapted  Solution 


The  final  section  presents  the  culmination  of  this  effort,  with  unsteady  solutions 
of  the  turbine  flow  field  made  using  periodic  grid  adaptation.  Starting  with  the 
coarse  grid  described  previously,  solutions  were  obtained  with  a  single  level  of  grid 
adaptation.  The  1:1  blade  count  ratio  was  again  used  to  reduce  the  computational 
cost  while  serving  to  highlight  the  improvement  in  the  solution  with  the  refined  grid 
solution.  The  adapted  solution  was  started  from  a  previously  converged  50  iteration 
per  passing  turbulent  solution  and  allowed  to  re-converge  to  a  periodic  solution.  As 
described  in  the  preceding  section,  a  study  was  performed  to  determine  the  number 
of  iterations  required  for  time  step  independence;  all  solutions  were  obtained  with 
two  orders-of-magnitude  drop  in  the  residuals.  Gradients  of  both  turbulent  eddy- 
viscosity  and  relative  Mach  number  were  used  to  trigger  grid  adaptation,  which 
occurred  20  times  per  blade  passing. 

Figures  6.46  -  6.49  show  the  changes  in  the  computational  grid  and  turbulent 
eddy  viscosity  at  four  times  during  a  blade  passing  period  for  the  adapted  solution. 
The  absolute  Mach  number  and  static  pressure  contours  for  two  of  these  times  are 
shown  in  Figs.  6.50  and  6.51.  Changes  in  both  the  vane  and  rotor  flow  fields  are 
quite  evident  as  the  adapted  portion  of  the  grid  tracks  the  location  of  the  vane 
wake.  As  the  rotor  passes,  the  accelerating  flow  around  the  suction  surface  thins, 
and  eventually  cuts,  the  vane  wake.  The  downstream  portion  of  the  wake  then  forms 
a  vortex  which  is  deposited  on  the  rotor  pressure  surface  and  remains  coherent  as  it 
is  convected  through  the  rotor.  It  is  clear  that  the  adaptation  routine  has  located 
both  the  region  of  rapid  acceleration  along  the  uncovered  portion  of  the  vane,  as 
well  as  the  pressure  surface  vortex. 

Figures  6.52  and  6.53  show  the  impact  of  the  grid  refinement  on  the  number 
of  iterations  required  to  obtain  a  time  step  independent  solution.  While  the  coarse 
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grid  solution  readied  an  asymptotic  solution  with  only  200  iterations  per  passing, 
the  adapted  solution  requires  800  iterations  per  passing.  This  result  is  expected 
after  the  experience  gained  through  computing  the  circular  cylinder  shedding  with 
the  finer  grid  in  an  earlier  section. 

Comparing  run  times  between  the  coarse  grid  and  adapted  solution,  there  is 
an  approximately  four-fold  increase  in  run-time  for  a  single  level  of  grid  adaptation. 
Coarse  grid  solutions  require  approximately  6  minutes  per  period  for  the  200  iteration 
per  period  case,  and  nearly  80  minutes  per  period  for  the  adapted  solution  with  800 
iterations  per  period.  The  number  of  points  in  the  solution  went  from  approximately 
5500  in  the  coarse  grid,  to  17,000  in  the  adapted  solution.  As  a  point  of  comparison, 
a  uniformly  fine  grid  would  contain  nearly  26,000  nodes  and  require  an  estimated 
116  minutes  per  period  for  the  800  iteration  per  passing  solution.  A  summary  of  the 
adapted  grid  solution  results  is  presented  in  Table  6.4. 


Adapted  Solutions 

real  time 

pseudo-time 

average 

CFL 

CPU  time 

iterations 

iterations 

iterations 

Number 

per  period 

per  period 

per  period 

per  time  step 

max  avg  min 

(min) 

100 

810 

8.1 

176  16.2  0.46 

30.2 

200 

1050 

5.3 

88  7.1  0.23 

29.6 

400 

1610 

4.0 

44  3.7  0.12 

47.0 

800 

2612 

3.3 

22  1.9  0.06 

76.4 

Table  6.4:  Results  of  iteration  count  per  period  tests;  adapted  grid,  2  orders-of- 
magnitude  residual  drop. 
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Figure  6.53:  Variation  in  unsteady  drag  coefficient  on  adapted  solution  of  turbine 
stage  with  different  iteration  counts. 
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Chapter  7 


CONCLUSIONS 


A  solution-adaptive  scheme  for  solving  the  quasi-three-dimensional  Navier-Stokes 
equations  through  turbomachinery  cascades  has  been  described.  The  method  uses 
central  differencing  of  both  the  inviscid  and  viscous  terms  and  a  one-equation  turbu¬ 
lence  model  for  the  Reynolds  stresses.  Time  integration  is  performed  in  an  implicit 
fashion  through  the  use  of  a  dual  time  stepping  approach.  The  computational  mesh 
is  composed  of  micro-blocks  of  quadrilateral  elements,  arranged  in  a  purely  unstruc¬ 
tured  fashion.  The  mesh  is  adapted  periodically  as  the  solution  evolves,  using  local 
cell  division  and  recombination  to  enhance  grid  resolution  at  user-selected  features 
in  the  flow  field.  This  chapter  summarizes  the  work  and  describes  the  conclusions 
that  have  been  derived  from  it.  Finally,  a  number  of  recommendations  are  presented 
for  improvements  to  the  code  which  might  be  explored  in  the  future. 

7.1  Summary 

The  aim  of  this  effort  was  to  produce  a  computational  tool  to  permit  the  rapid  com¬ 
putation  of  unsteady  flows  in  modern  transonic  turbomachinery.  The  end  product 
has  shown  the  viability  of  using  solution-adaptive,  unstructured  grid  methods  for 
unsteady  flow  computations.  Considerable  reduction  in  computer  time  is  possible 
through  the  use  of  local  cell  division  and  recombination  to  distribute  cells  only  where 
they  can  be  the  most  useful.  Computer  time  has  been  lowered  further  through  the 
implementation  of  the  dual  time  stepping  approach,  which  allows  an  explicit  time- 
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marching  method  to  be  used  as  a  driver  for  a  fully  implicit  scheme.  This  has  allowed 
the  use  of  convergence  acceleration  techniques,  while  maintaining  a  temporally  accu¬ 
rate  calculation. 


7.2  Conclusions 

The  following  section  summarizes  some  of  the  major  findings  of  this  effort: 

Dual  Time  Stepping  Technique.  The  flow  equations  are  recast  in  a  manner 
which  has  allowed  an  explicit  multi-stage  Runge-Kutta  scheme  to  be  used  as  the 
driver  for  a  fully  implicit  time  integration  method.  The  standard  time  derivative 
term  in  the  flow  equations  is  replaced  with  a  second-order,  backwards  difference  in 
time,  which  is  treated  as  an  additional  source  term.  A  pseudo-time  variable  is  then 
introduced  as  a  new  integration  variable  in  the  Runge-Kutta  scheme.  Because  the 
equations  are  no  longer  integrated  in  physical  time,  convergence  techniques  may  be 
used;  local  time  stepping  and  implicit  residual  smoothing  have  been  implemented. 
The  implicit  formulation  presented  here  is  slightly  different  from  the  methods  of 
Jameson  [32]  and  Weiss  and  Smith  [76]  found  in  the  literature,  with  a  fully  implicit 
treatment  of  the  vector  of  dependent  variables.  This  new  scheme  has  a  lower  opera¬ 
tion  count  than  either  method  and  does  not  require  the  limitation  on  the  pseudo-time 
step  size  required  by  Jameson.  A  significant  reduction  in  the  number  of  iterations 
and  computational  time  required  for  a  converged  periodic  solution  has  been  demon¬ 
strated  as  compared  to  a  solution  obtained  with  a  purely  explicit  method. 

Tests  have  been  performed  to  evaluate  the  performance  of  the  method  in  terms 
of  the  effect  of  the  magnitude  of  the  residual  drop  between  increments  in  physical 
time  and  the  effect  of  the  number  of  real  time  steps  per  period.  In  these  tests,  the 
results  have  been  shown  to  be  relatively  insensitive  to  the  residual  drop,  with  a  large 
change  in  computer  time  found  between  the  various  cases.  Results  are  considerably 
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more  sensitive  to  the  selection  of  the  number  of  iterations  used  per  period.  For 
most  solutions,  a  threshold  value  of  iteration  count  exists  above  which  only  small 
improvements  in  the  solution  appear  as  the  real  time  step  size  is  reduced.  Below  this 
threshold,  even  the  gross  unsteady  flow  features  are  poorly  predicted.  The  required 
iteration  count  per  period  is  also  shown  to  be  strongly  affected  by  the  relative  grid 
density.  Both  globally  refined  and  adapted  grids  demonstrate  a  requirement  for  an 
increased  number  of  iterations  per  period  to  obtain  a  temporally  converged  solution. 

One-Equation  Turbulence  Model.  The  effects  of  turbulence  have  been 
modeled  in  the  current  effort  with  the  one-equation  model  of  Spalart  and  Allmaras 
[71].  Since  the  model  contains  both  convective  and  diffusive  terms,  the  full  unsteady 
evolution  of  the  turbulent  eddy  viscosity  may  be  calculated  and  its  effect  on  the 
blade  row  interaction  evaluated.  This  model  has  been  implemented  with  an  integra¬ 
tion  scheme  similar  to  that  used  for  the  flow  equations  and  has  proven  to  be  quite 
robust.  By  including  the  turbulent  eddy  viscosity  in  the  viscous  time  step  limit  used 
to  compute  the  local  pseudo-time  step,  source  terms  in  the  model  may  be  evaluated 
explicitly  and  sub-iterations  of  the  turbulence  model  are  generally  not  required.  The 
turbulence  model  has  been  shown  to  provide  quite  reasonable  agreement  with  theo¬ 
retical  turbulent  boundary  layer  velocity  profiles  and  skin  friction  distributions,  even 
for  meshes  with  few  points  in  the  laminar  sublayer.  Compressor  cascade  computa¬ 
tions  show  the  prediction  of  premature  separation  on  the  blade  suction  surface  and 
highlight  the  difficulty  the  model  exhibits  for  cases  with  a  large  adverse  pressure 
gradient.  Blade  row  interactions  within  a  turbine  stage  have  been  computed  both 
with  and  without  the  turbulence  model,  as  well  as  for  the  case  where  the  turbulent 
eddy  viscosity  is  not  passed  between  the  blade  rows.  A  comparison  of  these  cases 
has  shown  that  the  propagation  of  eddy  viscosity  generated  by  the  upstream  blade 
row  has  only  a  small  effect  on  the  computed  lift  and  drag  coefficient  distributions. 
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Semi-Structured  Grids.  The  semi-structured  grid  approach  uses  micro- 
blocks  of  structured  cells  arrayed  in  a  globally  unstructured  fashion.  The  major 
benefits  of  this  grid  structure  are  the  large  reduction  in  the  memory  required  to 
store  the  grid  connectivity  array  and  the  lower  computational  time  resulting  from 
fewer  indirect  memory  references.  A  primary  reason  for  the  selection  of  this  type  of 
grid  structure  was  the  speed  at  which  grid  adaptation  could  be  performed.  Because 
the  logic  surrounding  the  adaptation  must  only  be  applied  to  the  various  patches, 
rather  than  to  each  individual  cell,  the  time  required  to  divide  or  fuse  a  patch  is 
greatly  reduced;  grid  adaptation  typically  occurs  in  less  than  the  time  of  a  single 
pseudo-time  iteration.  Finally,  the  semi-structured  approach  should  allow  a  natural 
implementation  of  the  multigrid  and  parallel  processing  techniques  described  in  the 
following  section. 

Blended  O-H  grids.  The  use  of  the  blended  0-H  grids  has  proven  quite 
successful  in  the  computation  of  turbomachinery  cascades,  for  both  compressors  and 
turbines.  Many  of  the  problems  that  occur  when  gridding  blunt  leading  and  trailing 
edges  have  been  removed,  while  maintaining  the  flexibility  of  an  H-grid  structure 
in  the  main  portion  of  the  domain.  Incorporation  of  the  multi-block  grid  structure 
has  further  improved  the  grid  through  the  removal  of  much  of  the  grid  skew  that 
occurs  at  the  point  where  more  than  two  0-  and  H-grid  blocks  meet.  Difficulties 
still  exist  in  the  generation  of  initial  grids  with  extremely  sheared  grids,  particularly 
for  turbine  cascades  with  their  high  turning  angles. 

Feature  Detection.  The  adaptation  of  the  computational  grid  is  performed 
to  locate  various  features  within  the  solution  flow  field.  Detection  of  these  features 
is  provided  in  two  ways.  A  strong  feature  detector  is  used  both  to  locate  the  regions 
of  large  pressure  gradient  found  near  shocks,  and  to  avoid  generation  of  additional 
grid  refinement  within  the  laminar  sublayer  of  a  turbulent  velocity  profile.  The 
smooth  feature  detector,  which  finds  regions  of  smooth,  but  rapid,  changes  in  the  flow 
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field,  uses  a  statistical  approach  based  on  undivided  differences  of  one  or  two  user- 
selected  flow  parameters.  Typically,  Mach  number  and  density  are  used  as  detection 
parameters,  but  turbine  blade  row  interaction  computations  have  shown  that  the 
turbulent  eddy  viscosity  provides  a  superior  method  to  track  the  propagation  of  the 
wake  through  the  downstream  blade  row. 

Smoothed  Interpolation.  When  a  micro-block  of  grid  cells  is  divided,  one 
method  to  locate  the  newly  generated  grid  points  is  by  simply  averaging  the  locations 
of  the  adjacent  points.  An  alternate  method  described  here  calculates  the  node 
locations  based  on  the  grid  stretching  currently  existing  with  the  patch.  This  method 
of  cell  division  has  proven  quite  effective  in  lowering  the  overall  level  of  grid  stretching 
error  in  the  mesh.  Reducing  the  stretching  error  is  of  particular  importance  in  the 
highly  stretched  regions  of  the  grid  near  blade  surface,  where  accurate  calculation 
of  the  viscous  terms  is  critical.  Maintaining  the  parent  patch’s  grid  line  distribution 
in  the  divided  patches  allows  the  stretching  ratio  to  be  halved  each  time  a  patch  is 
divided,  further  improving  the  accuracy  of  the  solution. 

7.3  Recommendations  for  Future  Work 

No  CFD  code  can  really  ever  reach  the  state  where  one  might  say  that  all  the 
developmental  work  is  truly  “finished.”  While  technology  improves  in  all  areas, 
progress  in  computational  algorithms,  turbulence  models,  and  computing  paradigms 
continues  to  drive  the  evolution  of  CFD  tools.  A  number  of  these  issues  are  discussed 
below: 

Differencing  of  Convective  Fluxes.  The  current  implementation  of  the 
convective  flux  differencing  uses  a  central  differencing  approach.  While  this  method 
has  a  number  of  advantages,  its  primary  disadvantage  is  the  degree  to  which  shocks 
are  smeared,  requiring  second-  and  fourth-order  damping.  Monotonic  methods, 
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which  difference  the  fluxes  along  the  characteristic  directions  in  the  flow,  have  a 
superior  ability  to  resolve  shocks,  normally  within  one  or  two  grid  points.  These 
methods  could  be  implemented  by  forming  the  additional  terms  from  the  monotonic 
scheme  as  a  replacement  for  the  existing  damping  terms,  while  retaining  the  current 
central  differencing  of  the  inviscid  fluxes. 

Turbulence  Model  Improvements.  The  Spalart-Allmaras  turbulence 
model  has  a  number  of  advantages,  both  in  its  treatment  of  the  turbulence  physics 
and  in  its  implementation.  While  this  has  made  it  attractive  for  this  effort,  it  has 
been  used  primarily  out  of  convenience.  Little  research  has  been  conducted  to  date 
to  construct  turbulence  models  that  are  inexpensive  to  use  and  simple  to  implement 
across  a  wide  range  of  unsteady  flow  problems.  Most  of  the  one-  and  two-equation 
models  that  are  available  for  engineering  use  do  well  in  low-speed,  steady  flows 
along  flat  plates  in  neutral  pressure  gradients.  Real-world  cases  involve  unsteady, 
transonic  flows  with  a  high  degree  of  wall  curvature  and  strong  pressure  gradients. 
It  should  be  noted,  however,  that  development  of  more  advanced  models  is  also 
driven  by  the  availability  of  detailed  test  data  in  these  more  complex  flow  regimes; 
data  acquisition  and  model  development  will  advance  hand-in-hand. 

Extension  to  Three-Dimensions.  While  the  use  of  the  quasi-3D  equations 
includes  much  of  the  flow  physics  found  in  modern  turbomachinery  designs,  a  number 
of  effects  are  lost  without  a  full  3-D  implementation.  These  3-D  effects  include  shock 
sweep,  endwall  boundary  layers,  tip  leakage,  and  other  secondary  flows.  Extension 
of  the  semi-structured  method  to  three  dimensions  should  be  fairly  straight-forward 
and  has  already  been  accomplished,  for  steady  flows,  by  Davis  and  Dannenhoffer 
[22]. 

Parallel  Processing  Implementation.  The  extension  of  the  current  work 
to  a  parallel  processing  environment  would  seem  to  be  a  natural  progression.  One 
of  the  primary  difficulties  of  porting  a  computer  code  to  either  a  fine-  or  coarse- 
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grained  parallel  processing  environment  is  the  proper  partitioning  of  the  solution 
domain  for  minimum  communication  overhead.  The  proper  goal  of  any  partitioning 
algorithm  is  to  minimize  the  length  of  the  boundary  at  the  interface  between  parti¬ 
tions.  The  current  implementation  has  an  advantage  over  traditional  unstructured 
codes  in  that  it  has  already  been  partially  partitioned  into  the  individual  patches. 
Partitioning  the  entire  flow  field  would  be  done  by  combining  all  the  patches  within 
some  level  of  parent  patch.  This  should  give  the  lowest  partition  boundary  length 
because  of  the  high  perimeter-to-area  ratio  of  the  logically  square  patch  shape.  Once 
domain  decomposition  is  accomplished,  each  partition’s  grid  and  flow  variables  are 
distributed  to  the  various  compute  nodes,  which  then  compute  independently.  Com¬ 
munication  of  the  boundary  data  could  occur  at  the  end  of  each  Runge-Kutta  stage, 
or  only  after  some  fraction  of  the  stages,  with  the  boundary  data  frozen  otherwise. 

Multigrid.  A  significant  advantage  of  the  dual  time  stepping  method  is  the 
ability  to  use  techniques  to  accelerate  solution  convergence.  Another  method  that 
might  be  implemented  is  multigrid,  where  a  series  of  successively  coarser  grids  is 
used  to  rapidly  propagate  changes  in  the  convective  fluxes  throughout  the  domain. 
In  the  present  semi-structured  implementation,  these  coarse  grids  may  be  formed 
in  one  of  two  ways.  The  first  method  is  to  remove  alternate  grid  lines  within  a 
patch  for  those  patches  at  the  finest  level.  Additional  levels  of  grid  coarseness  may 
be  obtained  using  the  same  method  applied  within  the  parents  of  divided  patches. 
Multigrid  is  particularly  useful  for  removing  low  frequency  errors  from  the  solution 
and  therefore  may  become  less  effective  as  the  size  of  the  time  step  in  physical  time 
is  reduced. 
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Appendix  A 


DERIVATION  OF  QUASI-3D  EQUATIONS 


This  appendix  presents  the  derivation  of  the  quasi-3D  form  of  the  Navier-Stokes 
equations.  As  described  in  Chapter  2,  these  equations  represent  a  two-dimensional 
formulation  along  a  blade-to-blade  streamsurface.  The  three-dimensional  effects  of 
streamsheet  thickness  variation  and  radius  change  are  then  added  and  appear  as 
source  terms  in  both  the  convective  and  diffusive  terms  of  the  equations. 

Repeating  the  conservation  equations  given  in  three-dimensional  form  in 
Eqns.  2.1  -  2.5  : 

4-  /  /  [  pdV  +  £  pVr  ■  ndS  =  0,  (A.l) 

dtJ  J  JV  JSy 

-riff  pVrdV  +  £  pV(Vr-n)dS  =  £  (— p  +  t)-MS ,  (A.2) 

dt  J  J  JV  JSy  J  sV 

J  JvEdV  +  jf  E(Vr-n)dS  =  <2  -  j£  (tV  -  pV)  ■  MS.  (A.3) 

In  order  to  expand  the  terms  in  these  equations  in  the  (m,  6 )  coordinate  system 
shown  in  Fig.  2.2,  vector  identities  are  required  for  this  system.  Using  the  method 
of  Anderson  [3,  pp.  193-197],  with: 

X\  =  r,  hi  =  h, 

x2  =  0,  h2  =  r, 

x3  =  m,  h3  =  1, 
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the  following  identities  may  be  obtained  (dropping  the  Ar  terms  since  vr  —  0  in  all 


1  d<j>„  06  A  1  06. 

v  <p  =  T~orlr  +  o  lm  ^ - 

h  Or  Om  r  00 


v ' a ^ k\'L{rhAm) + ~ls{hAi)  • 


where  6  is  an  arbitrary  scalar  and  A  is  an  arbitrary  vector.  The  velocities  V  and  Vr 
are  expressed  as: 


T  —  Vrnpm  T 

Id-  —  Vml'm  +  'UJO'l'Oi 


'here  w$  =  v$  —  rO  . 

Expressing  the  continuity  equation  in  the  equivalent  vector  differential  form: 


£^(V.i>)=!  +  V.(,V>)  =  0. 


Substituting  A  =  pVr  into  Eqn.  A. 5,  and  bringing  r  and  h  inside  the  temporal 
derivative  (since  neither  are  functions  of  time): 


S. :{rhp )  +  ~{rhpvm)  +  ^-(hpwe)  =  0. 


Similarly,  the  momentum  equation  may  be  expressed  as: 


DV  _  -  = 

',w  =  -v?+v'T' 


A. 10) 
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Expanding  the  LHS  of  this  equation  and  considering  only  the  m  and  9  components 
yields: 


DV  dV  -  — 

=  p-  +  ,K.W 
dV  fvodv„ 

'  ve  dwe 


,  „  -  ,  dvm  rm  \  „ 

f’-dl  +  p{7  -df  +  v^-vmV,r' 


, -  ,  dwe  rm\  „ 

+  p  7¥  +  ”-S  +  "'”"7 


(A.ll) 


where  the  following  notation  is  defined  [15]: 


1  dr 
r  dm  ’ 


1  dh 

h  dm 


(A. 12) 


Adding  zero  to  both  components  of  the  equation  in  the  form  of  the  continuity  equa¬ 
tion  and  combining  terms  yields: 


DV 
'  Dt 


+ 


Jt{rhpVm)  +  -L{rhp<)  +  Uhf,VM) ~  T )  Th 

jt(r2hpv,)  +  d-(r2hpvmv,)  +  ^(rhpweve)  -  pvevjD  7. 


(A.13) 


The  RHS  of  Eqn.  A. 10  is  also  expanded  using  the  method  outlined  in  Anderson  [3, 
pg.  196].  The  components  of  the  total  stress  tensor  II  =  —pt>ij  +  t  are  given  by: 


nr 

uB 


Ur 


P  T  e$$  &mrn)  •> 
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P  T  (2e#$  C ram  ^rr)  i 
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P  T  ^ rr  , 


(A.14) 
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nr0  n#r  — 

n  ffm  ffm(?  —  f-l &0ra  5 
n  mr  —  nrm  —  /J,6mr. 


The  expressions  for  the  strain  components  in  the  stress  tensor  are: 


m  h  ’ 

1  dvg  rm 

~~aa  +  vm  5 

r  dO  r 

dvm 

dm  ’ 

h  dr  r  h 
1  dVrn  d  V6_ 

r  dO  dm  r 


'A. 15) 


The  components  of  the  total  stress  tensor  are  then  given  by: 


2  [  h 
-P  +  -fi  2vm - 


U  _  ( 1  dye 
h  l  r  dO 


,2  / 1  dve  rm 

-p  +  r  2  Uw  +  ''”“ 


m  \  dvm 

r  J  dm  1 

dvm 

h  dm  5 
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( 1  dvQ 


0/3  Vm 
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(A.16) 
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Finally,  taking  V  •  II  and  rearranging  yields: 


v  •  n  =  —  Vp  +  v  •  t 

=  (  Tmm ))  +  -ggirTme)  +  (p  —  Tee)—  +  (p  — 


+  (; t<-r2hT"'e)  +  le{rh(-p+T‘ 


ie_ 
rh ' 


\  h"m  J 

'rT’~h  J  Hi 

(A. 17) 


where 


0  hm  2  n 

2 /ivTO—  -  — AtfcJ, 


=  2p- 


^1  <9ng 
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8m  -  I"®’ 


TmQ 


(dve  1 

^am  +  r  de 


'  m  \ 

-  ve — ), 

r 


(A.18) 


and  the  dilatation  is  given  by: 


( dvm  ldve 

e~  {l^  +  7~de  +v" 


V  h  N 
7 ~ +  X; 


(A.19) 


Combining  Eqns.  A. 13  and  A. 17  gives  the  two  momentum  equations: 

~(rhpvm)  +  -^{rh{pv2m  +  p  -  rmm))  +  h(pvmwe  -  TmB))  = 

=  r/i  ^(pu|  +  p  -  Tge)7^-  +  (p  -  ^  (A. 20) 

—(rhpvg)  +  -^{rh(pvmve  -  Tmg))  + —(h(pweve  +  p- Tee)) - 

=  -rh(pvevm  -  rme)  — .  (A. 21) 

r 
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By  multiplying  both  sides  of  the  0-direction  momentum  equation  by  r,  and  using 
the  chain  rule,  the  source  term  on  the  RHS  of  the  equation  may  be  absorbed  to  give: 


3  3  3 

— ( r2hpve )  +  —(r2h(pvmve  -  rmg))  +  —(rh(pweve  +  p 


=  0.  (A. 22) 


The  energy  equation  may  be  expressed  in  vector  differential  form  as: 


+  V  •  (EVr)  =  -V  •  q  +  V  •  (n  •  V). 


(A. 23) 


Using  the  same  technique  as  used  in  the  continuity  equation,  the  LHS  is  expanded 
to  yield: 


dE  -  1  ( dE  d  d  \ 

W  +  v ' (£K)  =  7h  br  +  1 [ThEVm)  +  do{ hEw,)  ■ 


(A. 24) 


Using  Fourier’s  Law,  q  —  —  kVX1,  and  neglecting  temperature  gradients  in  the  radial 
direction,  Eqns.  A.4  and  A. 5  give: 


-  _  1  /  d  ,  ,  dT,  d  /7  1  dr  A 

v'5=s  • 


(A. 25) 


Assuming  that  the  radial  shear  stresses  IIrm  and  Y[ro  are  zero  yields: 


U-V  =  A 


—  (^m(  P  T  ~mm  )  T  ^O^raB )  1>n 

+  {vmTm8  +  Vg(—  p  +  Tgg))  ig. 


(A. 26) 


(A. 27) 
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Combining  all  terms  and  simplifying,  the  final  form  of  the  energy  equation  is 
obtained: 


UrhE)  +  L{rh{E+p)vr' 

d 


*0) 


+  -qq  ( h(E  +  p)we  -  S4  +  rftp))  =  0, 


(A. 28) 


where 


D  9T 
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(A. 29) 
(A. 30) 


Finally,  the  continuity,  momentum,  and  energy  equations  may  be  expressed  in 
their  final  vector  form: 


dU_  dF_  8G_dR  dS 
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where 
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where  the  source  term  is: 


H2  —  (pVg  +  p  —  Tee) 

r 


(A.33) 


Appendix  B 


CHARACTERISTIC  BOUNDARY  CONDITIONS 

The  method  for  specifying  the  inflow  and  outflow  boundary  conditions  is  based  on  the 
development  of  non-reflecting  conditions  for  the  Euler  equations  described  by  Giles 
[27,  28].  By  using  characteristic  variables  to  specify  the  dependent  variable  values 
at  the  boundaries,  the  pressure  waves  generated  by  the  initial  solution  transients  are 
allowed  to  exit  the  domain  without  reflection.  This  results  in  faster  convergence,  an 
ability  to  place  computational  boundaries  closer  to  the  solid  bodies,  and  a  greater 
tolerance  of  the  flow  solver  to  poorly  specified  initial  conditions.  This  appendix 
discusses  the  theory  behind  these  boundary  conditions  and  describes  some  of  the 
special  considerations  required  by  turbomachinery  problems. 

B.l  Review  of  Characteristic  Theory 

Development  of  characteristic  theory  is  based  on  a  consideration  of  the  Euler  equa¬ 
tions  written  in  the  conservation  form  described  in  the  previous  Appendix.  However, 
the  development  in  this  Appendix  will  neglect  the  change  in  radius  and  streamsheet 
thickness.  The  assumption  is  made  that  the  inflow  and  outflow  boundaries  occur  in 
regions  of  nearly  constant  radius  and  streamsheet  thickness,  reasonable  assumptions 
for  axial  turbomachinery,  but  invalid  for  centrifugal  machines. 

The  Euler  equations  may  be  transformed  into  an  equivalent  form  based  on  the 
primitive  variables  p,vm,vg,p.  These  equations  are  linearized  by  considering  only 
small  perturbations  and  neglecting  the  higher  order  terms,  yielding  the  following 
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The  characteristic  variables  ci,  c2,  c3,  c4  represent  an  entropy  wave,  a  vorticity  wave, 
and  downward-  and  upward-running  pressure  waves  respectively;  c  is  the  isentropic 
sound  speed.  The  ci,C2,C3,c4  values  are  the  amplitude  of  these  waves,  and  p ,  um, 
he,  p  are  the  perturbations  from  the  uniform  flow  used  for  the  linearization. 

This  relation  is  used  to  develop  boundary  conditions  by  implementing  the 
conditions  as  modifications  to  the  flux  variables  at  the  boundaries,  rather  than  as 
a  specification  of  the  actual  boundary  values.  For  an  inflow  boundary,  either  three 
or  four  of  the  characteristic  variables  must  be  specified,  depending  upon  the  inflow 
Mach  number.  At  the  outflow  boundary,  the  specification  of  one  characteristic  may 
be  required,  again  depending  upon  the  Mach  number.  Characteristic  variables  that 
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are  not  specified  are  extrapolated  from  within  the  flow  domain:  the  mechanism 
which  allows  pressure  waves  to  exit  the  domain  at  either  boundary. 

If  the  change  in  the  primitive  variables  is  known  at  some  point,  then  the  change 
in  the  characteristic  variables  may  be  found  using  the  a  relation  similar  to  Eqn.  B.2: 
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with  the  inverse  given  by: 
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Here,  bc\,  bc2,  bc3,  bc4  are  the  calculated  characteristic  variable  amplitudes  and  bp, 
bvm,  bve,  bp  are  the  actual  changes  in  the  primitive  variables  from  which  the  changes 
in  the  dependent  variables  may  be  derived. 

The  basic  technique  in  formulating  these  boundary  conditions  is  to  first  choose 
an  appropriate  residual  to  be  reduced,  then  find  the  resultant  changes  in  the  char¬ 
acteristic  variables  required  to  achieve  this  reduction.  Changes  to  the  dependent 
variable  values  at  the  boundaries  are  then  calculated  using  the  relations  described 
above.  These  residuals  are  based  on  the  user-specified  flow  conditions  for  the  par¬ 
ticular  case  to  be  computed  and  must  be  chosen  with  regard  to  the  type  of  flow 
(subsonic  or  supersonic,  inflow  or  outflow)  that  is  expected  at  the  boundary.  The 
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following  sections  describe  these  various  choices  and  their  respective  implementa¬ 
tions. 

B.2  Inlet  Boundary  Conditions 

At  the  inflow  boundary,  at  least  three  of  the  four  characteristics  are  incoming,  with 
the  direction  of  the  fourth  determined  by  the  inflow  axial  Mach  number.  In  the 
special  case  of  a  supersonic  inlet  relative  flow  with  a  subsonic  axial  Mach  number, 
which  is  typical  of  transonic  compressors,  additional  constraints  are  placed  on  the 
boundary  condition  by  the  interaction  between  the  inlet  flow  and  the  geometry  of 
the  cascade.  The  following  section  discusses  the  implementation  of  these  various 
possibilities. 

B.2.1  Subsonic  Inlet  Flow 

For  a  subsonic  inlet  relative  flow,  the  inlet  flow  is  purely  subsonic  and  a  fairly 

standard  boundary  condition  may  be  implemented.  In  this  case  the  residuals  to 

be  driven  to  zero  are  based  on  the  inlet  entropy,  inlet  flow  angle,  and  stagnation 
enthalpy: 

Ri  =  p(S-S oo), 

R2  =  pc(ve  -  tan(aoo)um), 

R3  =  p(H-H, oo),  (B.5) 

where  S  is  an  entropy  function  given  by: 

S  =  log(7 p)  -  7  log(p),  (B.6) 
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and  H  is  the  stagnation  enthalpy: 


H  =  E  +  p.  (B.7) 

Nondimensionalizing  by  the  inlet  velocity  and  static  density,  the  specified  values  of 
these  functions  become: 
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Following  the  method  of  Giles  [28],  a  single  iteration  of  a  Newton- Ralphson 
procedure  is  performed  to  determine  the  changes  required  in  the  characteristic  vari¬ 
ables  to  drive  these  residuals  to  zero: 
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where  the  following  expressions  may  be  obtained: 
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1  0 


0  1  —  2  tan(aco) 

i  M,  |(1  +  Mm) 


(B.9) 


In  the  above  relation,  Mm  and  Mg  are  the  Mach  numbers  in  the  m  and  6  direc¬ 
tions  respectively.  Inverting  this  matrix,  the  required  changes  in  the  characteristic 
variables  are  obtained: 


1  +  Mm  +  Mg  tan(o:c 


1  +  Mm  +  Mg  tan(aco)  0  0 

tan(aco)  1  +  Mm  tan(a0 
— ir  -2  Mg  2 


R2  ■  (B.10) 


The  change  in  the  fourth  (outgoing)  characteristic  variable  is  determined  by 
extrapolating  its  value  from  a  point  adjacent  to  the  boundary,  using  the  values 
computed  during  the  inviscid  flux  integration: 


Sc4  =  -pc8vm  +  dp, 


(B.H) 


where 


=  (SUi-v^UJ/p, 

=  {SUs-vgSU^/p, 


=  (7  —  1)  (8U4  —  vm8\Ji  —  vgSUs  +  — (1^  ■  (B.12) 
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Here  the  SU  values  are  evaluated  at  the  next  point  away  from  the  inlet  boundary. 
The  changes  in  the  primitive  variables  may  be  obtained  by  applying  Eqn.  B.4: 


6p  = 

Sv  — 
8ve  = 

8p  — 


c 2  L 
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—8c\  +  —(8cz  T  8c4 ) 


2pJS°3  SC 
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pc  ’ 


~(8cs  +  8c4). 


(B.13) 


Finally,  from  these  values,  the  required  changes  in  the  dependent  variables  at 
the  inlet  boundary  are  found: 


SUi  =  8p , 

SU2  =  p8vm  -J-  vm8p , 

8U3  =  pSvg  +  VgSp, 

8U4  =  -^-r  +  vM  +  veSUs-hvl  +  vDSlh.  (B.14) 

7-1  2 

The  radius  and  streamsheet  contraction  terms  must  also  be  included  in  these 
dependent  variable  changes;  it  is  also  important  to  remove  these  terms  when  extrap¬ 
olating  8c4  in  Eqn.  B.ll. 


B.2.2  Supersonic  Inlet  Flow 

For  a  supersonic  inlet  flow  condition,  with  the  axial  Mach  number  still  subsonic,  an 
alteration  of  the  inlet  boundary  condition  is  required.  In  a  transonic  compressor 
cascade,  it  is  the  interaction  between  the  inlet  flow  and  the  geometry  of  the  cascade 
which  sets  the  inlet  flow  angle;  the  mass  flow  is  set  by  the  exit  static  pressure.  This 
is  referred  to  as  the  unique  incidence  condition  [18].  Because  of  this  condition,  the 
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inlet  flow  angle  can  no  longer  be  specified;  instead,  the  inlet  relative  Mach  number 
or  the  tangential  velocity  (wheel  speed)  must  be  given.  The  inlet  flow  angle  is 
then  determined  after  the  cascade  flow  has  converged.  This  condition  is  also  the 
reason  why  a  transonic  fan  has  such  a  small  variation  in  mass  flow  over  its  operating 
characteristic;  there  is  only  a  small  range  of  acceptable  inlet  flow  angles  between 
the  point  at  which  the  cascade  chokes  and  where  the  bow  shock  spills,  causing  the 
cascade  to  stall.  Many  researchers,  including  Giles,  miss  this  point  and  incorrectly 
specify  the  inlet  flow  angle.  Use  of  this  condition  may  yield  a  convergent  solution, 
however  it  does  deny  some  of  the  important  physical  phenomena  at  work. 

The  specification  of  the  supersonic  inlet  boundary  condition  proceeds  in  a 
fashion  similar  to  that  developed  for  the  subsonic  condition.  In  this  case,  the  resid¬ 
uals  to  be  driven  to  zero  are  based  on  the  inlet  entropy,  tangential  velocity,  and 
stagnation  enthalpy: 


R-i  =  p(S  -  Soo), 

R2  =  pc{ye  -  v&oa), 

Rs  =  p{H-H oo).  (B.15) 


The  procedure  outlined  previously  is  used,  yielding  the  following  differences  from 
the  subsonic  results: 
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Inverting  this  matrix  gives  the  required  change  in  the  characteristic  variables: 
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The  boundary  condition  specification  is  completed  using  a  procedure  identical  to 
that  used  for  the  subsonic  condition. 


B.2.3  Pure  Supersonic  Inflow 

For  an  inlet  where  the  flow  is  purely  supersonic,  i.e.  has  a  supersonic  axial  Mach 
number,  all  four  characteristic  variables  must  be  specified.  Assuming  that  this  is  also 
a  steady  flow  specification,  there  will  be  no  change  in  the  characteristic  variables  from 
the  values  used  to  originally  specify  the  initial  conditions  of  the  problem.  Therefore, 
the  inlet  condition  is  quite  simple: 

8cx  =  8c2  =  8c3  =  8ca  =  0,  (B.18) 


or  equivalently: 


8b\  =  8U2  =  8U3  =  8U4  =  0. 


(B.19) 


This  is  the  boundary  condition  formulation  used  in  the  supersonic  flat  plate 
boundary  layer  cases  described  previously. 
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B.3  Outflow  Conditions 


At  the  outflow  boundary,  at  least  three  of  the  four  characteristics  are  outgoing,  with 
the  direction  of  the  fourth  determined  by  the  exit  axial  Mach  number.  The  following 
sections  describe  these  two  possibilities: 

B.3.1  Supersonic  Outflow? 

If  the  exit  axial  Mach  number  is  supersonic,  then  the  fourth  characteristic  is  also 
outgoing,  and  no  specification  of  the  exit  quantities  is  allowed.  The  changes  in  the 
flow  variables  at  the  exit  boundary  are  found  by  a  simple  extrapolation  from  the 
fluxes  calculated  at  the  adjacent  grid  point: 


=  6Ux  .  , 

Ww 

=  SU2 

SUZex.t 

=  su3 

exit 

SU4exit 

=  SU4  .  . 

exit~ 

(B.20) 

The  subscript  exit  indicates  points  along  the  exit  boundary,  and  exit~  indicates 
points  just  in  from  the  boundary.  Again,  it  should  be  noted  that  the  ratios  of  the 
streamsheet  contraction  and  radius  must  be  included  in  this  evaluation. 

B.3. 2  Subsonic  Outflow 

For  a  subsonic  exit  Mach  number  condition,  the  downstream  static  pressure  is  typ¬ 
ically  specified.  This  is  equivalent  to  the  setting  the  cascade  pressure  ratio,  often 
prescribed  in  compressor  and  turbine  cascade  calculations.  In  this  case,  the  fourth 
characteristic  variable  is  incoming  and  must  therefore  be  determined.  Here,  the 
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residual  to  be  driven  to  zero  is  based  on  the  difference  between  the  actual  and  spec¬ 
ified  static  pressure.  The  fourth  characteristic  variable  is  given  by: 

6c4  -  -2 (p  -  pexit),  (B.21) 

where  pexu  is  the  user  specified  value. 

The  other  three  characteristic  variable  values  are  determined  using  the  simple 
extrapolation  described  above  to  give: 

<5  ci  .  =  Sc.\  , 

Sc2ex„  =  8c2cx.t_, 

Sc3  ,  =  6c3  .  .  (B.22) 

Having  the  values  for  the  change  in  all  four  characteristics  at  the  exit  boundary,  the 
changes  in  the  dependent  variables  may  be  calculated  as  before. 

B.3.3  Mixed  Exit  Condition 

A  mixed  subsonic-supersonic  exit  condition  has  also  been  implemented  for  use  in 
the  supersonic  boundary  layer  calculations  used  for  flow  solver  validation.  In  this 
case,  a  search  out  from  the  wall  along  the  exit  boundary  is  performed  to  locate  the 
first  point  at  which  the  flow  becomes  supersonic.  The  supersonic  points  are  treated 
with  the  extrapolation  condition,  while  subsonic  points  use  the  pressure  at  the  first 
supersonic  point  to  set  their  exit  pressure. 
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Appendix  C 


INTERACTIVE  GRID  GENERATOR 


This  appendix  describes  the  grid  generator  used  to  provide  the  initial  computational 
grids  for  cascade  calculations.  As  discussed  in  Chapter  5,  the  grids  produced  are  in  a 
multi-block  format  and  include  both  single  and  multiple  blade  row  cases.  Addition¬ 
ally,  the  grid  generator  is  responsible  for  the  creation  of  the  streamsheet  thickness 
and  radius  variation  data  files.  The  gridding  tool  uses  a  graphical  user  interface 
(GUI)  which  allows  the  user  to  access  virtually  all  of  the  parameters  needed  to 
generate  a  multi-block  grid.  The  following  sections  describe  the  geometrical  data 
requirements,  the  multi-block  grid  format,  the  elliptic  and  algebraic  grid  generation 
algorithms,  and  the  graphical  interface.  This  appendix  serves  as  a  brief  introduction 
to  the  grid  generator  and  is  not  intended  as  an  user’s  guide,  this  will  be  left  to  a 
later  publication. 


C.l  Multi-Block  Grid  Structure 

The  grids  constructed  for  this  thesis  use  a  multi-block  O-H  grid  structure  decompo¬ 
sition  of  the  computational  domain.  Using  the  combination  of  these  two  grid  types 
allows  the  user  to  take  advantage  of  the  relative  advantages  of  each  type.  The  topo¬ 
logical  decomposition  of  a  typical  grid  is  shown  in  Fig.  C.l.  An  O-grid  is  wrapped 
around  the  blade  to  use  its  superior  ability  to  produce  grid  lines  which  are  nearly 
normal  to  the  blade  surface.  Further,  O-grids  allow  leading  and  trailing  edges  with 
small  radii  to  be  gridded  with  minimal  grid  skew.  The  O-grid  is  split  into  two  halves, 
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Figure  C.l:  Logical  topology  for  multi-block  grids. 
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divided  along  the  periodic  boundaries  which  extend  from  the  leading  and  trailing 
edges.  These  two  blade  grids  are  blended  into  a  passage  grid  which  is  composed  of 
multiple  H-grid  blocks.  These  blocks  consist  of  a  main  passage  grid  (labeled  A  in 
Fig.  C.l)  and  optional  blocks  upstream  (blocks  B  and  D)  and  downstream  (blocks 
C  and  E)  of  the  two  O-grid  halves.  During  the  grid  generation  procedure  described 
below,  the  user  has  control  over  the  number  of  grid  lines  in  each  of  these  blocks. 
The  H-grid  structure  was  selected  for  the  passage  grid  to  facilitate  the  construction 
of  the  inlet,  exit,  and  sliding  interface  boundary  conditions,  which  occur  along  the 
right  and  left  edges  of  the  H-grid  block.  Periodic  boundaries  in  the  H-grid  portion 
occur  between  the  bottom  of  block  B  and  top  of  block  D,  and  similarly  between 
blocks  C  and  E.  All  H-grid  blocks  are  stored  together  in  a  global  H-grid  domain 
with  counters  used  to  locate  the  corners  of  the  various  blocks.  These  counters  are 
written  to  a  separate  file  used  by  the  flow  solver  to  reconstruct  the  H-grid  blocks 
and  position  the  O-grid  halves. 

C.2  Geometry  Definition  Requirements 

In  order  to  facilitate  the  integration  of  the  computational  tool  described  in  this 
thesis  into  the  design  system  used  in  the  US  Air  Force  Wright  Laboratory,  Aero 
Propulsion  and  Power  Directorate,  a  standard  data  format  was  adopted.  This  format 
was  taken  from  the  geometry  definition  files  generated  by  the  in-house  fan  and 
compressor  design  program  UD0300.  This  program  uses  a  streamline  curvature 
method  to  develop  blade  shapes  from  user  input  performance  parameters;  a  limited 
description  of  this  program  is  available  in  [40,  77]. 

Two  types  of  geometrical  data  files  are  used  by  the  grid  generator  in  the  cre¬ 
ation  of  a  quasi-3D  blade-to-blade  grid  surface.  The  first  is  a  definition  of  the  blade 
shape  on  a  series  of  streamsurfaces  ordered  from  hub  to  tip.  These  shapes  are  true 
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blade  surface 


Figure  C.2:  Surface  data  used  to  construct  blade  trailing  edge. 

3-D  surfaces  which  use  a  local  coordinate  system  centered  at  the  stacking  axis  of 
the  blade.  Typical  UD0300  program  output  provides  80  points  on  each  of  the  blade 
surfaces,  with  another  30  points  used  to  describe  the  leading  edge  radius;  the  trailing 
edge  is  assumed  to  be  blunt.  The  grid  generator  modifies  the  blunt  trailing  edge 
condition  by  fitting  a  circular  arc  within  the  wedge  angle  formed  by  the  last  surface 
segments  on  each  of  the  two  sides  of  the  blade,  as  shown  in  Fig.  C.2. 

A  second  data  file  is  used  to  position  the  blade  shapes,  each  having  their  own 
local  coordinate  system,  within  the  global  coordinate  system  of  a  multi-blade  row 
machine.  This  file  is  simply  the  meridional  ( r  —  z )  plane  axisymmetric  computational 
mesh  which  evolves  in  the  streamline  curvature  algorithm  of  the  UD0300  program;  an 
example  mesh  for  a  transonic  compressor  is  shown  in  Fig.  C.3.  This  mesh  describes 
a  series  of  concentric,  annular  streamtubes  through  the  machine  and  is  arranged 
such  that  the  mass  flow  is  equal  in  each  of  the  streamtubes  and  remains  constant 
through  the  length  of  the  machine.  This  results  in  a  zero  flow  condition  across  the 
streamtube  boundaries,  which  is  one  of  the  fundamental  assumptions  of  the  quasi- 
3D  formulation.  This  global  coordinate  system  file  also  allows  multiple  blade  row 
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Figure  C.3:  Meridional  mesh  for  a  transonic  compressor. 

cases  to  be  generated  by  providing  the  necessary  positioning  information  for  each  of 
the  separate  blade  surface  files. 

The  streamtube  thickness  and  radius  variation  information  is  generated  by 
combining  these  two  file  types.  For  a  grid  developed  along  a  particular  streamsurface, 
the  radius  variation  is  taken  from  the  radial  values  along  the  streamsurface;  thickness 
values  are  determined  by  taking  the  difference  in  radius  between  the  two  adjacent 
streamsheets. 


C.3  Grid  Generation 

The  generation  of  the  computational  grid  proceeds  in  two  phases.  The  first  is  the 
definition  of  the  computational  boundaries  and  the  various  grid  parameter,  followed 
by  a  series  of  smoothing  steps  to  achieve  an  acceptable  final  grid. 

The  computational  domain  is  defined  by  the  blade  surfaces,  the  inlet  and  exit 
boundaries,  and  the  periodic  boundaries  up-  and  downstream  of  the  leading  and 


191 


trailing  edges  respectively.  The  user  has  control  of  the  length  and  angle  of  the 
up-  and  downstream  sections,  as  well  as  the  boundary  between  the  H-  and  O-grid 
portions  of  the  grid.  Counter  widgets  on  the  graphical  interface  are  used  to  set  the 
number  of  grid  lines  in  each  direction  of  the  various  grid  blocks.  These  counters 
are  operated  by  repeated  mouse  clicks  to  increment  or  decrement  to  the  desired 
value,  or  by  typing  in  the  number  directly.  The  current  value  of  the  patch  size  is 
used  when  incrementing  or  decrementing  counter  values  to  insure  that  connectivity 
requirements  are  satisfied  between  the  various  grid  blocks.  Finally,  slider-bar  widgets 
are  provided  to  set  the  coefficient  of  the  exponential  stretching  functions  applied 
along  and  normal  to  the  blade  surface,  on  the  periodic  boundaries,  and  across  the 
passage  between  blades.  The  grid  is  updated  with  these  new  values  in  real-time, 
allowing  the  user  to  immediately  see  the  impact  of  any  changes  in  these  variables. 

Once  the  grid  definition  parameters  are  selected,  an  initial  grid  is  constructed 
by  simply  connecting  the  appropriate  points  along  the  various  boundaries.  Because 
of  the  complexities  of  the  geometry  in  a  typical  turbomachinery  cascade,  this  grid 
is  unacceptable  for  computation  and  must  be  smoothed.  Two  methods  are  used 
to  smooth  this  initial  grid,  providing  a  grid  with  smoothly  varying  cell  areas  and 
minimal  cell  skewing.  The  first  is  the  elliptic  method  of  Sorenson  [70],  which  may 
be  applied  throughout  the  entire  H-grid  portion  of  the  grid  or  only  in  the  center 
passage  H-grid  section.  A  simple  Laplacian  smoother  is  then  available  to  uniformly 
smooth  the  entire  grid.  Each  of  these  smoothing  algorithms  is  available  individually, 
or  may  be  used  in  a  combined  mode  with  an  iteration  of  each  performed  alternately. 
A  typical  grid  smoothing  session  proceeds  with  a  near  convergence  of  the  elliptic 
smoother  followed  by  a  number  of  Laplacian  smoothing  iterations  to  equalize  the 
cell  sizes  along  the  H-  and  O-grid  boundaries. 
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C.4  Graphical  Interface 


The  gridding  tool  graphical  interface  was  designed  to  allow  the  user  to  access  vir¬ 
tually  all  of  the  parameters  needed  to  generate  the  grid.  Each  of  these  quantities  is 
available  and  may  be  changed  at  any  point  in  the  gridding  procedure;  any  parameter 
change  is  immediately  incorporated  in  the  user’s  view  of  the  grid.  An  example  of 
the  GUI  is  shown  in  Fig.  C.4. 

The  GUI  is  divided  functionally  into  five  separate  windows,  or  panels:  the 
large  grid  window  at  the  bottom  left  displays  the  current  view  of  the  grid  during 
generation;  the  main  control  panel  along  the  top  is  used  to  control  the  file  I/O  and 
contains  the  dials  which  are  used  to  zoom  and  translate  the  grid  in  the  grid  window. 
This  panel  also  contains  buttons  which  determine  the  portions  of  the  grid  that  are 
visible,  and  controls  the  application  of  the  smoothing  operators.  The  panel  at  the 
upper  right  is  used  to  set  the  main  geometrical  parameters  of  the  computational 
domain,  such  as  the  up-  and  downstream  periodic  boundary  lengths  and  angles  and 
the  spacing  between  blade  rows.  Counters  for  the  row  number,  number  of  passages  in 
a  particular  blade  row,  and  the  patch  size  are  also  located  here.  The  last  two  panels 
along  the  right  side  contain  the  counters  and  sliders  for  the  grid  dimensions  and 
boundary  stretching  values  for  the  H-  and  O-grid  portions  of  the  grid  respectively. 
Since  the  edges  of  all  blocks  in  the  grid  must  contain  an  integer  number  of  patches, 
the  dimension  counters  are  incremented  or  decremented  by  the  current  patch  size. 

The  GUI  has  been  implemented  using  the  widget  set  from  the  “panel  library” 
of  the  NASA-developed  FAST  post-processor  [75].  Unfortunately,  this  library  limits 
use  of  the  grid  generator  to  Silicon  Graphics  workstations;  conversion  of  the  interface 
to  the  Motif  widget  set  is  planned  for  the  future. 
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Figure  C.4:  Grid  generator  graphical  user  interface;  2:3  blade  ratio  turbine  case  shown. 


C.5  Example  Mesh  Generation 

In  order  to  illustrate  the  use  of  the  program,  an  example  grid  generation  session  is 
presented.  This  case  is  a  two  blade  row  transonic  turbine  stage  with  a  1:1  blade 
count  ratio. 

After  entering  the  names  of  the  meridional  mesh  and  two  blade  streamline 
coordinate  files,  the  user  sees  the  “grid  box,”  which  consists  of  outlines  of  the  blade 
surfaces  and  the  edges  of  the  computational  domain.  If  the  user  is  generating  an  0-H 
grid,  the  outer  edge  of  the  O-grid  surrounding  the  blade  surface  is  also  displayed. 
As  mentioned  previously,  sliders  in  the  GUI  allow  the  user  to  adjust  the  upstream 
and  downstream  grid  boundary  angles  and  lengths,  the  distance  between  blade  rows, 
and  the  size  of  the  O-grid  oval  around  the  blades.  The  grid  window  for  this  initial 
grid  box  view  is  shown  in  Fig.  C.5. 

The  next  step  is  to  fill  the  grid  boundaries  with  the  unsmoothed  initial  grid.  At 
this  point,  the  user  selects  the  number  of  grid  points  in  each  of  the  blocks  of  a  multi- 
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Figure  C.6:  Initial  grid  in  “filled”  state. 


block  grid,  and  the  stretching  functions  applied  along  the  various  grid  boundaries. 
Grid  points  are  located  along  the  boundaries  in  accordance  with  the  selected  count 
and  stretching  values  and  the  connectivity  requirements  between  the  various  grid 
blocks.  This  is  termed  the  “filled  state”  and  is  shown  in  Fig.  C.6. 

This  initial  grid  is  now  smoothed  with  the  application  of  the  elliptic  and  Lapla- 
cian  smoothing  functions.  Typically,  elliptic  smoothing  is  first  performed  to  achieve 
a  nearly  converged  grid  solution  in  the  H-grid  blocks.  The  Laplacian  function  is 
then  applied  to  blend  the  interface  between  the  O-  and  H-grid  portions  of  the  grid. 
The  GUI  contains  buttons  for  each  of  these  smoothers;  each  button-push  causes  one 
iteration  of  the  selected  smoothing  algorithm  to  be  performed.  A  close-up  of  the 
region  between  blade  rows  after  the  elliptic  smoothing  is  shown  in  Fig.  C.7;  the  final 
grid,  after  application  of  the  Laplacian  smoother,  appears  in  Fig.  C.8. 

The  final  step  in  the  process  is  the  selection  of  an  output  filename  and  grid  file 
type.  The  grid  data  is  written  in  the  format  used  by  the  Plot3D  and  FAST  programs 
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Figure  C.7:  Grid  after  elliptic  smoothing. 


Figure  C.8:  Final  grid  after  Laplacian  smoothing. 
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[12,  74,  75]  and  may  be  written  in  either  a  formatted  or  unformatted  form.  The  user 
is  also  given  the  option  to  write  a  small  file  containing  all  the  parameters  used  to 
generate  the  grid  in  its  pre-smoothed  state.  This  file  can  be  saved  and  read  in  later 
to  restart  the  grid  generation  procedure.  The  grid  output  routine  also  writes  the 
files  which  contain  the  radius  and  streamtube  thickness  variations,  the  blade  surface 
spline  data,  and  the  multi-block  connectivity  information. 
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